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ABSTRACT 


The  design  of  a  composite  panel  requires  some  way  of  finding  the 
minimum  thickness  laminate  which  will  withstand  the  load  requirements 
without  failure.  The  mathematical  complexity  of  this  problem  dictates 
the  use  of  non-linear  optimization  techniques.  Although  there  are 
sophisticated  optimization  programs  available  capable  of  solving  for  the 
ply  ratios,  these  programs  are  not  often  used  in  preliminary  design 
because  they  require  a  large  computer  and  some  knowledge  of  the 
program's  operation.  As  an  alternative,  specialized  laminate 
optimization  programs  were  developed  which  are  compact  and  efficient 
enough  to  run  on  microcomputers.  Only  stresses  at  a  point  and  inplane 
loads  and  deflections  are  considered.  The  programs  are  simple  to  use 
and  require  no  knowledge  of  optimization.  Techniques  are  developed  in 
this  thesis  that  find  minimum  thickness  laminates  with  either  ply  ratios 
or  ply  angles  as  design  variables.  In  addition,  a  method  is  presented 
for  finding  the  optimimum  orientation  for  the  axis  of  symmetry  of  an 
orthotropic  laminate.  The  orthotropic  laminate  program  uses  an 
approximate  failure  theory,  as  suggested  by  Tsai,  that  has  been  found  to 
speed  computations  dramatically.  £ — """ 

Many  test  cases  were  run  with  these  programs  to  demonstrate  the 
weight  savings  possible  over  quasi-isotropic  laminates.  Of  particular 
interest  is  performance  of  the  laminates  under  multiple  independent 
loads.  Initial  orientations  for  the  programs  to  operate  on  were 
studied,  and  0/90/45/-45  laminates  were  found  to  be  an  effective 
starting  point  for  design. 

The  approximate  failure  criterion  made  analytic  investigations  of 
optimized  laminates  possible.  A  method  of  plotting  maximum  strain 


energy  density  as  a  function  of  the  shear-stress-free  laminate 
orientation  is  derived  t  demonstrate  how  the  laminates  adapt  to 
multiple  design  load  requirements  in  the  optimization  process.  Also,  an 
optimality  criterion  is  derived  which  is  satisfied  by  each  ply  group  at 
the  minimum  thickness  condition. 


I.  INTRODUCTION 


Background 

Almost  any  introductory  text  on  composite  materials  will  make  a 
statement  to  the  effect  that  one  of  the  principle  advantages  of 
composites  is  the  possibility  of  tailoring  the  laminate  to  suit  the 
structural  requirements.  By  using  the  directional  nature  of  the 
material  to  advantage,  highly  efficient  structures  should  be  possible. 
Yet,  except  for  uniaxial  loads,  no  suggestion  is  made  for  selecting 
these  tailored  laminates.  The  omission  is  not  accidental,  but  is  due  to 
the  difficulity  of  converting  the  equations  for  laminate  analysis  into 
equations  for  laminate  design. 

When  sizing  an  isotropic  plate,  the  orientations  and  the  number  of 
plies  at  each  orientation  can  be  variable.  Although  analysis  equations 
for  finding  the  response  of  a  given  laminate  are  well  known,  these 
equations  cannot  be  solved  to  yield  the  best  laminate  for  a  given  set  of 
requirements.  Besides  being  non-linear,  structural  design  requirements, 
such  as  strength,  are  stated  as  inequalities.  There  is  no  way  to  know 
how  to  assign  equalities  to  the  equations  and  solve  for  the  design 
variables.  We  cannot  tell  which  combination  of  requirements  will  be 
"critical"  for  the  best  design. 

A  common  approach  to  sizing  laminates  is  to  asume  the  plies  are 
acting  independently.  For  strength  requirements,  this  is  referred  to  as 
netting  analysis.  Although  in  general  this  is  a  bad  approximation, 
reasonable  results  can  be  obtained  for  0/90  laminates  with  no  shear. 

With  any  other  case,  such  as  additional  ply  groups,  off-axis  loads,  or 
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multiple  independent  loads,  netting  analysis  cannot  provide  a  reasonable 
basis  for  design.  The  plies  within  a  laminate  interact  in  a  complex 
manner  and  cannot  be  considered  independent.  Because  of  the 
interaction,  there  are  no  simple  formulas  for  proper  sizing,  nor  is 
intuition  a  reliable  guide. 

Non-linear  optimization  techniques  developed  over  the  last  20  years 
provide  a  sound  mathematical  basis  for  laminate  sizing.  The  techniques 
should  not  be  thought  of  as  the  final  step  in  design,  used  to  shave  off 
a  couple  percent  of  weight,  but  as  the  starting  point  of  design. 
Optimization  can  be  applied  to  any  design  constraint  that  can  be 
mathematically  modelled.  Constraints  may  include  stiffness,  strength, 
stability,  minimum  gage,  and  dynamic  response.  In  this  thesis,  the 
author  has  chosen  to  work  only  with  strength  constraints.  Besides  being 
an  essential  element  of  design,  it  is  one  of  the  few  constraints  which 
can  be  described  as  a  point  problem,  assuming  loads  do  not  change  as  the 
laminate  changes.  The  optimal  laminate  will  be  considered  as  one  with 
minimum  thickness,  and  thus  weight.  For  constraints  such  as  stability, 
the  optimizer  must  be  coupled  to  a  structual  analysis  method,  such  as 
finite  elements,  in  order  to  describe  the  geometry  and  boundary 
condition  influences.  The  assumption  that  optimization  for  strength  can 
be  dealt  with  as  a  point  problem  is  completely  valid  only  for  a 
determinate  structure.  The  optimization  procedure  will  have  to  be 
coupled  to  a  structural  analysis  code  in  some  iterative  process  in  order 
to  properly  size  indeterminate  structural  elements.  Still,  the  simple 
methods  and  programs  presented  here  should  be  of  aid  in  much  of  the 
initial  design  process. 

The  role  of  optimization  is  particularly  inportant  when  designing 
for  multiple  loading  conditions.  A  wing  panel  must  sustain  several 


different  flight  conditions,  as  well  as  ground  loads.  Not  only  are  the 
magnitudes  of  these  loads  changing  with  time,  but  the  orientation  of  the 
load  principle  axes  may  also  change.  For  directional  materials,  it  will 
often  be  convenient  to  think  in  terms  of  shear-free  loads  and  an  angle 
that  transforms  the  loads  to  the  laminate  axis  system.  Because  of  the 
laminate's  anisotropic  strength,  changes  in  the  principle  axis  leads  to 
a  problem  that  does  not  exist  for  isotropic  materials}  it  is  impossible 
to  pick  a  severest  load  condition  by  inspection  and  size  the  laminate  to 
that  load  alone.  In  fact,  there  may  not  be  any  single  severest 
condition.  For  a  minimum  weight  laminate,  some  of  the  plies  may  be  near 
failure  for  one  load,  while  other  plies  are  critical  for  a  different 
load.  One  result  of  this  added  complication  is  that  optimization 
results  cannot  be  tabulated  in  a  design  manual.  There  is  no  way  to 
characterize  all  the  possible  loading  combinations  into  a  finite  set  of 
graphs.  Instead,  the  computer  must  become  an  integral  part  of 
preliminary  design. 

If  optimization  is  so  valuable  to  the  design  of  composite 
laminates,  why  isn't  it  in  common  useage?  After  all,  the  basic  methods 
of  non-linear  optimization  are  well  developed  and  can  handle  much  more 
complex  problems  than  sizing  a  laminate.  Indeed,  laminate  sizing  is  a 
comparatively  well  behaved  problem,  with  typically  only  a  few  design 
variables  and  constraints.  Part  of  the  answer  may  be  the  reluctance  to 
use  a  complex  and  general  code  requiring  a  main-frame  computer.  In 
addition,  there  may  be  some  lack  of  confidence  in  the  procedure.  This 
thesis  presents  some  specialized,  user-friendly  codes  which  can  be  run 
on  microcomputers  at  the  designer's  desk.  Hopefully,  by  having  a  desk¬ 
top  computer  that  only  requires  the  user  to  respond  to  some  simple 
prompts  for  input,  further  application  of  optimization  will  be 


encouraged. 

The  potential  for  applying  optimization  techniques  to  composites 
has  not  escaped  the  attention  of  other  authors.  At  least  2  programs 
exist  in  a  documented,  publicly  available  form.  One  is  COMAND  by 
Vanderplaats  [1]  which  couples  a  laminate  analysis  program  to  an 
existing  general  optimization  code,  also  by  the  same  author.  Maximum 
strain  failure  criteria  are  used,  and  minimum  values  of  the  A  matrix 
components  can  be  entered  to  account  for  minimum  stiffness  requirements. 
Another  program  was  written  by  Khot  [2].  Instead  of  a  direct  numerical 
optimization,  this  program  relies  on  the  assumption  that  strain  energy 
density  will  be  equal  for  all  ply  groups  as  the  laminate  approaches 
minimum  thickness  [3],  An  iterative  procedure  for  adjusting  the  number 
of  plies  is  derived  from  the  assumed  optimality  condition.  The  program 
also  includes  an  approximate  buckling  constraint,  based  on  "smeared" 
laminate  properties.  The  optimization  routines  are  coupled  to  a  finite 
element  code  to  update  the  stress  state  as  the  composite  panels  change. 
Neither  of  these  programs  meets  the  requirement  for  simplicity  of  use 
which  is  the  goal  of  this  thesis. 

Without  a  numerical  optimization  program,  the  minimum  thickness 
laminates  can  still  be  studied  if  there  is  only  one  free  variable,  such 
as  the  best  angle  in  an  angle-ply  laminate.  Some  of  these  one¬ 
dimensional  searches  are  presented  in  [4].  This  reference  is  notable 
because  it  includes  the  approximate,  strain-sphere  failure  criterion 
discussed  later  in  this  thesis. 

The  programs  written  in  the  course  of  this  work  are  all  in  BASIC. 
The  particular  computers  were  chosen  somewhat  arbitrarily,  but  the  codes 
should  be  readily  transferable  to  other  computers  with  a  minimum  of 
change.  Optimization  with  the  quadratic  failure  criteria  with  a 


complete  set  of  laminate  property  outputs  requires  about  12  kilo-bytes 
of  memory.  The  angle  optimization  can  be  attached  for  about  2k  more 
memory.  A  simplified  version  based  on  an  approximate  failure  criteria 
fits  in  less  than  6K.  Programs  have  been  written  for  the  Timex-Sinclair 
1000  [5],  the  Epson  HX-20  [6],  and  the  Texas  Instruments  CC-40  .  These 
last  2  microcomputers  were  picked  because  they  offer  true  desk-top 
capability;  the  original  goal  of  the  project. 


Laminate  Theory 


The  developement  of  the  laminate  plate  theory  equations  will  follow 
Tsai  and  Hahn  [7]  wherever  possible.  The  difference  will  be  that  vector 
notation  is  used  more  extensively  in  this  thesis.  The  plates  will  be 
subject  only  to  inplane  loads  and  deflections.  The  order  of  plies  in 
the  laminate,  or  stacking  sequence,  is  not  a  factor  in  the  optimization 
procedure.  However,  for  the  inplane  deflection  restriction  to  be  valid, 
the  actual  laminate  would  have  to  be  symmetric.  That  is,  for  any  ply  at 
orientation  6 ,  a  distance  Z  above  the  midplane,  there  is  a 
corresponding  ply  of  the  same  orientation  at  minus  Z.  For  these 
restrictions,  strain  is  a  constant  through  the  thickness  and  the  stress- 
strain  relation  is  simply 


N  *  [A]e 

where 


0) 


m 


(1), 

1=1  >  ''1 


Ajk  »  «£,  Q^'h 


(2) 


^-laminate  strain  vector 

ft-load  vector  in  terms  of  stress  resultants 

(1) 

Qj|<  -modulus  component  trarsformed  from  the  orientation 
of  the  i’th  ply  group 
m-nunber  of  ply  groups 
h.|  -thickness  of  the  i'th  ply  group 

Several  ways  exist  to  perform  the  transformations.  The  programs  listed 
in  Appendices  B-D  use  an  invariant  formulation  with  multiple-angle 
functions  as  given  in  reference  [7].  In  terms  of  engineering  constants, 
the  Q’s  are  given  by 


Qxx=mE> 


Qyy-mEy 


QXy=mvyEx 


Qss=  Es 


where 

m  =  0-  Vy)_1 

Ex  is  the  longitudinal  Young's  modulus,  vx  the  longitudinal  Poisson's 
ratio,  Ey  the  transverse  Young's  modulus,  and  vy  the  transverse 
Poisson's  ratio. 

The  axis  system  convention  is  shown  in  Figure  1.  x,  y,  and  s 
subscripts  denote  properties  in  the  ply  axis  system,  and  1,  2,  and  6 
denote  properties  in  the  laminate  axis  system. 

A  ply  group  will  be  defined  as  all  the  plies  of  a  particular 
orientation  and  material  (for  hybrids),  whether  or  not  they  are  actually 
adjacent  in  the  laminate.  In  the  optimization  procedure,  ply  group 
thickness  is  handled  as  a  continous  variable.  The  individual  ply  as  a 
discret  unit  is  ignored.  After  the  procedure  is  finished,  we  must 
divide  the  ply  group  thickness  by  the  thickness  of  an  individual  ply  and 
round-off  to  get  the  integer  mumber  of  plies  required.  A  logical  way  of 
rounding-off  must  be  a  topic  of  future  research.  For  now,  rounding-up 
should  be  assumed  for  all  ply  groups.  The  term  "ply  ratio"  will  also  be 
used.  This  is  the  ratio  of  a  particular  ply  group  thickness  to  the 
total  laminate  thickness. 

For  the  graphs  and  tables  presented  in  this  thesis,  the 
conventional  lamination  code  becomes  awkward.  Instead,  the  notation 

(0/90/t 45) 

refers  to  the  class  of  laminates  with  those  orientations,  with  ply  group 
thickness  determined  by  the  optimization  procedure.  Also, 

(0i  /901  /i45-j  ) 

refers  to  a  laminate  with  the  stated  orientations  and  equal  ply  ratios. 
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Failure  Criteria 


One  of  several  failure  criteria  could  be  selected  for  incorporation 
in  the  optimization  procedure  [8],  The  quadratic  tensor  polynomial  or 
Tsai-Wu  criterion  was  selected  because  it  fits  experimental  data  well 
[7]  and  because  it  reduces  the  number  of  constraints  as  compared  to 
maximum  stress  or  strain  criteria.  The  quadratic  failure  criterion  is 
based  on  fitting  an  ellipse  to  the  experimental  failure  strengths  of  a 
unidirectional  lamina.  The  form  of  the  equation  accounts  for 
interaction  between  the  stresses  causing  failure.  As  in  most  laminate 
failure  criteria,  each  ply  in  the  laminate  must  be  interrogated 
seperately  in  order  to  determine  if  failure  has  occured.  In  this 
thesis,  first-ply  failure  is  adopted,  in  contrast  to  a  progressive 
failure  model. 

The  quadratic  failure  criterion  takes  the  form 


+  F1°i  "  1  £  0  i.J  =  1.2,6 

The  F's  are  related  to  experimental  data  as  follows 

1 


xx  ~nr 


F  =1-1 

rx  x  X' 


1 


yy  yy  ' 


f  =1-1 

\y  y  r 


Fxy  =  Fxy  ^xx^yy 


where  X  -longitudinal  tensile  strength 


(4) 


(5) 


fc  --!■  —  >-  - 
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X'  -longitudinal  compression  strength 


Y  -transverse  tensile  strength 
Y'  -transverse  compression  strength 
S  -shear  strength 


F  -non-dimensional  interaction  term 
xy 


F^y  has  not  yet  been  accurately  measured  since  it  requires  a  reliable 
biaxial  stress  test.  From  geometric  bounds  and  by  analogy  to  isotropic 
materials  (Von  Mises  failure  theory)  a  value  of  -1/2  is  usually  taken, 
and  is  used  throughout  this  thesis. 

Stating  the  failure  criteria  in  terms  of  strain  is  convenient.  In 
strain  space  the  failure  envelopes  stay  fixed  even  if  the  ply  ratios  of 
the  laminate  are  changed.  The  strain  limits  of  a  ply  are  independent  of 
the  laminate  stiffness.  This  is  an  important  conceptual  simplification 
when  ply  ratios  are  variable.  The  failure  criterion  can  be  rewritten  as 


GjjVj  *  G(sf  '  0 


1,j=  1,2,6 


where  the  G's  are  found  by  applying  the  stress-strain  relations, 
assuming  linear  elasticity  to  failure.  Then 


GkJt  °  F1 

Gj  .  FjQ,j 


i » J  »k  ,2,  =  1,2,6 


The  G  and  F  matrices  can  be  transformed  for  off-axis  plies  by  a  second- 
order  tensor  transformation,  just  as  with  the  elasticity  components. 

The  linear  terms  of  the  equation  (G  vector)  are  transformed  by 

G^  =  P  +  q  cos  29 
G?  =  P  -  q  cos  20 

( 

Gg  =  q  sin  29 


P  =  1  {Gx  +  Gy)  :  q  =  \  (Gx  -  G  ) 

Figure  2  shows  the  failure  enveloped  for  a  0/90  laminate  of 
T300/5208  (Graphita/Epoxy) .  The  envelopes  are  actually  three- 
dimensional,  and  shear  strain  is  not  shown.  Only  the  region  enclosed  by 
both  ellipsoids  is  considered  safe. 

An  approximate  first-ply  failure  envelope  was  suggested  in 
reference  [7].  The  envelope  is  based  on  recognizing  there  is  a  first- 
ply  failure  domain  common  to  all  possible  ply  orientations,  and  thus 
independent  of  the  orientation  of  any  particular  ply.  Figure  3  shows 
failure  envelopes  for  several  orientations.  There  is  an  inner  envelope 
defined  by  the  0°and  90°plies,  within  which  no  failure  occurs  for  any 
possible  orientation.  Note  that  0°and  90°plies  do  not  always  define 
this  space  for  other  material  systems.  By  using  this  inner  envelope,  we 
have  a  failure  criterion  which  applies  to  the  laminate  as  a  whole,  and 
does  not  need  to  be  interrogated  on  a  ply-by-ply  basis.  It  is 
convenient  to  fit  some  analytic  surface  into  the  envelope.  Since 
tension  loads  are  of  primary  interest  in  this  work  (because  there  are  no 
stability  constraints),  a  sphere  centered  on  the  origin  was  selected  to 
give  a  conservative  approximation  of  the  inner  envelope.  The 
approximate  failure  criterion  can  then  be  written 


+  e2  +  1  zl  <  b2  (9) 

1  2  2  o  - 

The  sphere's  radius,  b,  can  be  set  equal  to  the  minimum  lamina  strain, 
ti. -an  directly  from  experimental  data.  The  criterion  will  be  referred 
to  as  a  maximum  strain-sphere. 

The  strain-sphere  criterion  will  not  be  acceptable  for  uniaxial 
laminates,  or  for  loads  in  the  3rd  quadrant  (compression-compression), 
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but  otherwise  is  of  some  value.  The  simplicity  of  the  criterion  more 
than  doubles  the  speed  of  the  optimization  algorithm.  For  optimization 
with  tension-tension  loads,  it  has  been  found  to  be  about  7% 
conservative,  as  compared  to  the  Tsai-Wu  criterion.  Thus,  for  quick 
answers,  the  approximation  is  adequate.  In  addition  to  allowing  for 
extra  fast  computation,  the  maximum  strain-sphere  is  simple  enough  to 
allow  analytic  investigations  of  the  optimization  process,  as  will  be 
discussed  in  later  sections. 
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Failure  Constraints  in  Design  Space 


Normally,  failure  envelopes  show  the  set  of  loads  (or  strains)  that 
can  be  sustained  by  a  particular  laminate.  For  design  purposes,  the  set 
of  laminates  that  can  sustain  particular  loads  would  be  more  desirable. 
Instead  of  stress  or  strain  coordinates  on  a  graph,  the  coordinates 
should  be  the  design  variables,  for  example,  ply  group  thicknesses. 
Unfortunately,  there  may  be  an  arbitrary  number  of  design  variables,  and 
therefore  dimensions  to  the  problem.  Therefore,  general  design  graphs 
cannot  actually  be  drawn,  but  the  concept  is  important  to  tinderstanding 
the  optimization  process. 

One  way  of  showing  the  set  of  laminates  that  could  sustain  a  given 
combination  of  loads  is  to  make  a  plot  which  divides  design  space  into 
two  regions;  a  region  where  the  laminates  would  not  fail  for  any  of  the 
given  loads  (called  the  feasible  region),  and  a  region  where  the 
laminates  would  fail  (called  the  infeasible  region).  Any  point  in 
design  space  defines  a  unique  laminate.  We  will  restrict  the  discussion 
to  taking  ply  group  thicknesses  as  the  only  design  variables.  The 
boundary  between  the  feasible  and  infeasible  regions  is  the  surface 
defined  by  the  the  failure  criteria  equations  when  made  into  an  equality 
and  plotted  as  functions  of  the  thicknesses.  With  the  quadratic  failure 
criterion,  we  can  write 


-T 

C(l) 


(10) 


where  the  subscript  L  designates  the  strains  associated  with  a 
particular  set  of  loads  and  superscript  P  denotes  a  transformation  from 
a  particular  angle.  Equation  (10)  can  be  shown  to  be  a  function  of  the 
h’s  (ply  group  thickness)  by  substituting 
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An  important  feature  of  working  in  design  space  is  that  the  constraint 
surfaces  for  more  than  one  set  of  loads  can  be  plotted  together.  The 
final  result  is  several  surfaces  in  the  design  space,  with  the  outermost 
surfaces  forming  the  boundary  between  feasible  and  infeasible  space 
(Figure  4). 

If  there  are  only  two  design  variables,  we  can  actually  draw  these 
design  graphs.  Figures  5,  6,  and  7  are  plots  of  the  constraint  curves 
for  a  0/90  laminate  under  a  single  biaxial  load.  The  three  figures  are 
for  three  different  failure  criteria.  To  define  the  feasible  region, 
the  maximum  strain  criterion  requires  the  number  of  surfaces  to  be  three 
times  the  number  of  ply  groups  times  the  number  of  independent  loads 
(only  4  curves  are  shown  in  Figure  5  because  shear  strain  is  zero  for 
the  particular  class  of  laminate  and  the  given  load).  The  quadratic 
criterion  requires  the  number  of  surfaces  to  be  equal  to  the  number  of 
ply  groups  times  the  number  of  loads.  Reducing  the  number  of 
constraints  speeds  the  optimization  procedure.  Speed  of  operation  is 
another  motivation  for  choosing  the  quadratic  criteria  for  the  majority 
of  additional  work.  The  approximate  strain-sphere  criterion  is  simplier 
yet,  with  only  a  single  surface  for  each  independent  load.  Because  it 
is  a  conservative  approximation,  only  limited  use  will  be  made  of  this 
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INFEASIBLE  SPACE 


FIGURE  5:  Failure  Constraints  for  Maximum  Strain  Criteria 
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II.  OPTIMIZATION  METHODS 


Ply  Ratios 


The  laminate  sizing  problem  can  be  stated  in  the  language  of 
optimization  theory  as  follows; 
find  min.  of  h 


where 


subject  to 


h  = 


ro 

Z 

i=l 


(12) 


where 


CP,L  «  0 

P,i  ■  1,2.. 

h.j  >  0 

L  =  1,2.... 

r  -+T 

,  r(P) ,  + 

CP.L  *  e(D 

lG  1  S(L) 

+ 


5(p)Te 


(L)  ‘ 


1 


ls<p>l 


S<p> 


-quadratic  failure. criteria  parameters 
transformed  from  the  orientation  of  ply 
group  P 

-linear  terms  of  failure  criteria  transformed 
from  the  orientation  of  ply  group  P 
-component  of  strain  due  to  loading  L 


Although  simply  stated,  there  is  no  simple  solution.  One  of  several  non¬ 
linear  optimization  methods  could  be  applied  to  the  problem.  A 
modification  of  the  method  of  feasible  directions  was  chosen  after 
examining  ways  to  speed  the  computations  enough  so  that  solution  on  a 
microcomputer  could  be  practical.  The  modification  of  the  method  makes 
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use  of  certain  closed  form  equations  at  intermediate  steps,  reducing  the 
number  of  calculations  needed.  The  algorithm  also  takes  advantage  of 
the  linearity  of  the  objective  function  in  terms  of  the  design 
variables.  This  simplification  also  speeds  up  the  algorithm  as  compared 
to  more  general  formulations. 

Although  nany  figures  in  this  section  show  the  optimization  process 
on  two-dimensional  graphs  in  design  space,  it's  important  to  realize 
that  sone  aspects  of  the  problem  may  not  be  evident  until  3  or  more 
dimensions  are  considered.  For  example,  the  constraints  may  form  long, 
narrow  valleys  that  the  search  method  must  follow  efficiently.  Because 
all  mathematics  are  derived  in  vector  form,  the  extension  to  higher 
dimensions  is  simply  a  matter  of  book  keeping  for  the  computer. 

Design  optimization  must  always  take  into  account  the  issue  of 
local  versus  global  minima.  From  optimization  theory,  if  the  feasible 
space  can  be  shown  to  be  convex,  then  there  is  only  a  global  minima  [9], 
An  informal  definition  of  convexity  is  that  any  two  points  in  the  space 
can  be  connected  by  a  straight  line  which  does  not  pass  out  of  the  space 
at  any  point.  The  intersection  of  convex  spaces  forms  a  convex  space 
[9],  Thus,  if  each  constraint  surface  is  convex,  then  there  is  only  one 
minima.  From  observation  of  actual  plots  for  cases  with  2  ply  groups, 
the  failure  constraints  of  composites  meet  this  requirement.  No  proof 
of  the  generality  of  this  observation  is  offered,  but  the  assumption 
that  the  optimization  leads  to  a  global  minima  from  any  starting  point 
will  be  accepted  in  this  thesis. 

Due  to  the  periodicity  of  trigonometric  functions,  there  will  not 
be  a  single  minima  when  angles  are  varied.  This  is  a  severe  handicap  to 
making  angles  a  design  variable. 


In  the  method  of  feasible  directions,  the  design  is  changed  so  that 
the  trajectory  in  design  space  follows  the  constraint  surfaces  along  a 
direction  that  decreases  the  objective  function  as  quickly  as  possible, 
but  never  leaves  the  feasible  region.  A  non-linear  constraint  cannot 
be  followed  continuously  because,  numerically,  the  algorithm  must  take 
finite,  linear  steps.  Therefore,  a  vector  is  found  which  both  decreases 
the  objective  function  and  does  not  violate  the  constraint  for  a  finite 
move.  The  trajectory  of  a  feasible  direction  algorithm  is  shown  in 
Figure  8. 

The  problem  with  this  method,  for  our  purposes,  is  that  finding  the 
distance  to  the  next  constraint  along  an  arbitray  vector  requires  a 
numerical,  one-dimensional  search.  Since  each  constraint  evaluation 
requires  forming  the  laminate  A  matrix,  inverting  the  matrix,  solving 
for  strains,  and  evaluating  the  failure  equations,  we  would  like  to 
reduce  the  number  of  iterations  required  for  this  search.  Some 
approximations  were  tried,  based  on  assuming  the  inverse  of  strain  to  be 
a  linear  function  of  ply  group  thickness.  These  were  meant  to  speed  the 
search,  but  were  not  found  to  be  completely  reliable.  Instead,  the 
method  was  modified  to  allow  for  larger  error  bands  in  the  numerical 
search. 

Briefly,  the  modification  consists  of  measuring  the  distance  across 
the  constraint  surface  "valley",  along  a  vector  on  which  the  objective 
function  is  a  constant.  This  restricts  the  method  to  problems  with  an 
objective  function  that  is  linear  in  terms  of  the  design  variables. 
Finding  this  distance  still  requires  a  numerical  one-dimensional  search, 
such  as  bisection,  but  now  the  error  band  can  be  quite  large,  reducing 
the  number  of  iterations  needed.  The  larger  error  band  is  allowable 
because  only  a  rough  measure  of  the  distance  across  is  needed,  whereas 
in  the  feasible  directions  method,  the  constraint  surface  must  be 


located  with  high  accuracy,  since  that  point  serves  as  the  starting 
coordinate  of  the  next  iteration  of  the  search.  We  assume  the  bottom  of 
the  "valley"  will  be  about  halfway  across.  From  the  halfway  point,  ply 
ratios  are  keep  constant,  and  the  total  thickness  of  the  laminate  is 
scaled  so  that  the  coordinates  in  design  space  rest  directly  on  the 
constraint  surface  defining  the  feasible  region.  The  scaling  operation 
is  based  on  recognizing  that  for  constant  ply  ratios,  strain  is 
proportional  to  total  thickness.  This  closed  form  equation  compensates 
for  the  error  band  of  the  numereial  search.  From  the  new  coordinate, 
the  procedure  repeats  until  changes  are  very  small,  or  a  new  search 
direction  cannot  be  found  (Kuhn-Tucker  conditions  for  optimality  [9]). 

A  possible  trajectory  for  the  modified  method  is  shown  in  Figure  9. 

The  constraint  that  thickness  be  greater-than  or  equal-to  zero  is 
known  as  a  "side  constraint".  These  linear  constraints  are  simple 
enough  to  be  handled  by  seperate  logic.  If  the  one-dimensional  search 
hits  a  side  constraint,  and  no  strength  constraints  are  violated  at  that 
point,  the  procedure  stops  on  the  h  =0  plane,  rescales  the  '  minate,  and 
proceeds  as  before.  Any  constraints  associated  with  a  zero  thickness 
ply  are  ignored.  Once  a  ply  is  set  to  zero  thickness,  it  is  never 
restored.  The  ability  to  completely  drop  a  ply  group's  constraints 
seems  to  be  unique  to  the  programs  developed  for  this  thesis. 

A  step-by-step  description  of  the  algorithm  will  be  presented, 
along  with  the  relevant  equations.  For  clarity,  the  variables  used  in 
this  section  will  not  always  be  identical  to  those  actually  used  in  the 
programs . 
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ajectory  for  Modified  Method 


1)  Laminate  Scaling 

Before  any  optimization  of  ply  ratios  can  be  considered,  we  must 
first  be  able  to  size  the  total  thickness  of  a  laminate  with  constant 
ply  ratios.  Strains  are  proportional  to  total  thickness.  This  is 
evident  by  writting  the  stress-  strain  relation  as 

e  =  £  [a*]  N  (13) 

where  J^a*jis  the  thickness  normalized  inverse  of  the  A  matrix. 

Instead  of  total  thickness,  it  is  more  convenient  to  use  the  change  in 
the  distance  from  the  origin  in  design  space  as  the  scaling  parameter. 
The  strain  proportionality  is  the  same  for  either  parameter  since 


h  .  mil-  A  .  -  UAh°)2  _ 

"  Zh*  "  *  r°  Eh°2 


(14) 


where  &  is  a.  proportional  change  of  the  individual  ply  thicknesses,  and 
r  is  the  distance  from  the  origin.  To  use  this  linear  relation,  a 
reference  strain  vector  is  calculated,  along  with  a  reference  r°  Then, 
as  long  as  ply  ratios  are  consta,  t,  strain  for  any  other  value  of  r  can 
be  found  from  the  equation 


e 


(15) 


where  the  superscript  o  refers  to  reference  conditions.  This  relation 
can  be  substituted  into  any  of  the  strain-space  failure  criteria.  With 
the  quadratic  criterion  we  have 

.(P) 


V  |G'p)|?.  ♦  f  G(P,V  =  1  -  e 

(L)  '  1  (l)  r  (L)  1 


(16) 
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where  e j  is  a  small  (10  )  offset  that  ensures  the  point  stays  slightly 


in  the  feasible  region  despite  any  numerical  error.  Solving  for  r 


where 


A 

B 

C 


=  1  -  e 


=  -G 


1 

(P>V  r° 


“*£>  |G 


(L) 

(P) 


£(l)  r 


The  value  of  r  should  be  calculated  for  every  possible  constraint.  The 
largest  resulting  value  corresponds  to  the  constraint  forming  the 
boundary  between  feasible  and  infeasible  space.  With  this  value  of  r, 
the  ply  group  thicknesses  are  scaled  according  to 


0 


where  again,  the  superscript  o  means  a  reference  condition. 


2)  Initial  Feasible  Point 


Thicknesses  are  first  set  to  a  large,  arbitrary  value,  to  be 
assured  of  starting  in  the  feasible  region.  The  program  sets  all  ply 
group  thicknesses  to  l/»^m  where  m  is  the  number  of  ply  groups.  Next, 
the  total  thickness  of  the  laminate  is  scaled  so  that  one  constraint  is 
critical  (Figure  10).  The  scaling  operate  is  given  above. 


3)  Active  Constraint  List 


At  any  step  in  the  optimization,  one  or  more  constraints  will  be 


active.  These  are  the  constraints  that  are  currently  near  critical  as 
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CONSTRAINT  SURFACES 


FIGURE  10:  Scale  Total  Laminate  Thickness 
Ply  ratios  constant 


defined  by 


eT  |G(P>(  ?  +  r(P)T"  , 

en  l  |b  I  e/n  +  G  x  -  1  >  e. 


where  a  value  of  0.05  has  been  found  to  work  well  for  e^.  Before 
finding  a  search  direction,  the  program  must  evaluate  this  equation  for 
all  values  of  P  and  L,  and  maintain  a  list  of  these  values  for  which  the 
constraint  is  active. 


4)  New  Direction 


We  need  to  find  a  vector  which  points  away  from  all  the  active 
constraints  and  is  parellel  to  the  constant  total  thickness  plane 
(Figure  11).  Components  of  the  gradient  vector  are  first  calculated  for 
each  active  constraint  according  to  the  equation 


-P,L  .  tT  .  (P)r  .  .(P)U 

)hi  (L)  ,G  ^(D.hi  G  e(L).hi 


where 


(L) ,hi  =  de2/dhi 
de6/dh. 


Since  the  applied  loads  are  independent  of  the  laminate  configuration, 
the  partials  of  strain  can  be  evaluated  from  the  stress-strain  relation 
as  follows; 


0  =  ^(1^) 

=  | A | ,hje  +  |A|e,hi 
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where 


I A  |  »hi  =  |Q°;| 

so  that 

t,h,  =  -| A*1 |  |Q(1)!e 


The  gradient  of  each  active  constraint  is  normalized  to  unit  length. 

The  individual  gradients  are  then  summed  and  the  result  is  normalized  to 
a  unit  length.  The  reason  for  summing  the  gradients  can  only  be 
visualized  in  3  dimensions.  Suppose  two  constraint  surfaces  meet  to 
form  a  valley,  and  the  objective  function  can  still  be  reduced  by 
following  the  valley  along  its  length.  If  only  one  constraint  were 
operated  on  at  a  time,  the  trajectory  would  bounce  inefficiently  back 
and  forth  between  the  surfaces.  By  summing  the  normalized  vectors,  a 
resulting  vector  that  points  down  the  valley  can  be  formed.  The 
negative  of  the  summation  will  point  into  feasible  space.  This 
resultant  vector  will  be  called  W. 

The  projection  onto  the  constant  thickness  plane  is  done  by  the 


double  cross-product 


z  =  n  x  (w  x  n) 


which,  by  a  vector  indenity  can  be  written 

z  =  w  -  (w  •  n/n 

where  n  is  the  unit  normal  to  the  plane  defined  by 


Z  h-  =  constant  (r 

i=l  1  v 

In  keeping  with  good  numerical  practice,  Z  is  also  normalized.  If  the 

length  of  Z  before  normalization  is  small  (  10  ^  )  then  W  and  n  must  be 

near  parallel.  This  would  indicate  that  a  minimum  has  been  reached  and 
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the  program  halts 


5)  Distance  to  Next  Constraint 


The  next  step  is  to  find  the  distance  along  Z  to  the  next 
constraint  (Figure  12).  A  bisection  method  is  used  for  the  one- 
dimensional  search.  The  vector  Z  describes  relative  changes  in  the  ply 
group  thicknesses.  Moving  a  scalar  distance  S  along  Z  changes  the 
thicknesses  according  to 

h  =  h°  +  SZ  (24) 

where  h°  is  the  vector  of  current  thicknesses  for  S=0.  Note  that  even 

though  the  individual  ply  groups  are  changing,  total  thickness  stays 
constant  along  Z.  The  program  will  need  to  be  able  to  quickly  calcuate 
the  A  matrix  as  ply  groups  change.  To  save  a  few  multiplications,  the 
programs  represents  A  as 

lAl  =  (A0I  +  s|Az|  (25) 

where 


*0, 


k=l 


m 

Azij  =k=l^ij  zk 


The  initial  bounds  on  the  bisection  search  are  S=0  and  S=Smax  where 
Smax  is  the  distance  to  the  nearest  h=  0  constraint.  Smax  is  calculated 
by  finding  the  largest  positive  value  of  the  equation 


^max  =  1  =  (26) 

The  usual  bisection  method  is  slightly  modified.  First,  instead  of 
trying  to  find  the  zero  of  a  single  equation,  we  must  evaluate  each 
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possible  constraint  to  find  the  boundary  between  feasible  and  infeasible 
spaces.  The  programs  in  the  appendices  contain  a  subroutine  which 
evaluates  the  constraints  and  returns  a  single  flag  with  the  value 
"FAIL"  if  a  single  constraint  is  violated,  and  "PASS"  if  no  constraint 
is  violated  (Cp)L<0  for  all  P,L)  The  second  feature  is  that  S=Smax  may 
be  in  the  feasible  region.  What  this  means  is  that  a  ply  group  can  be 
reduced  to  zero  thickness  without  violating  any  constraints.  If  this  is 
the  case,  the  program  updates  the  ply  group  thickness  vector  for  the 
point  S=Smax  and  rescales  the  laminate,  eliminating  constraints 
associated  with  the  zero  thickness  ply  group.  The  algorithm  then 
restarts  from  step  2.  If  S=Smax  is  not  feasible,  then  the  bisection 
continues  with  the  follow  steps: 

1)  Let  S1=0,  S2=Smax 

2)  Let  S=(Sl+S2)/2 

3)  Test  all  constraints  at  point  S 

4)  If  f lag="PASS"  then  S1=S 

If  flag="FAIL"  then  S2=S 

5)  If  S2-S1<10-^  then  search  direction  immediately  hits 
constraint.  This  indicates  the  minimum  has  been 
found. 

6)  If  (S2-Sl)/Sl>l/4  then  go  to  step  2.  Else  stop 
bisection  procedure 

Step  6  checks  to  see  if  the  error  with  which  the  distance  to  the 

constraint  is  known,  is  less  than  1/4  the  distance  across  the  "valley". 

The  1/4  is  arbitray,  but  gives  good  overall  convergence  of  the  algorithm 
with  a  minimum  number  of  bisection  iterations.  Note  that  for  each  value 
of  S  tested,  the  A  matrix  must  be  formed,  inverted,  strains  calculated, 
and  constraint  evaluated. 


5)  Rescale  Laminate 


Once  the  distance  to  the  next  constraint  is  known,  we  take  S=Sl/2. 
From  this  point  in  design  space,  the  total  thickness  is  reduced  by  the 
laminate  scaling  procedure  (Figure  13).  If  the  change  in  total 
thickness  is  small  (less  than  1/10  a  single  ply  thickness),  the 
algorithm  is  assumed  to  have  reached  a  minimum  and  halts.  If  not,  the 
algorithm  repeats  from  step  2.  The  loop  continues  until  one  of  the  halt 
conditions  is  reached. 

The  organization  of  the  program  is  shown  by  a  flowchart  in  Figure 
14.  The  flowchart  is  only  meant  to  be  an  aid  to  understanding  the  steps 
required.  The  interconnections  between  subroutines  in  the  actual 
programs  are  somewhat  more  complex. 

Table  1  gives  some  examples  of  the  convergence  rate  and  number  of 
inverse  A  matrix  evaluations  (the  most  time  consuming  step)  required  for 
the  optimization.  Times  are  given  for  a  ZX-81  computer  which  has  a  Z-80 
microprocessor.  An  iteration  is  counted  as  the  total  loop  from  step  2 
to  5.  Three  or  4  iterations  is  typical  unless  some  ply  groups  are  going 
to  zero  thickness,  which  counts  as  a  full  iteration. 


NEW  DESIGN 
COORDINATES 


FIGURE  1 3:  Rescale  Total  Thickness 
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Continuously  Variable  Angles 


The  most  obvious  approach  to  selecting  appropriate  ply  orientations 
is  to  let  the  computer  calculate  the  optimal  values.  There  is  no 
fundamental  reason  why  this  cannot  be  done,  but  there  are  some 
implementation  problems,  and  the  results  are  not  always  satisfactory. 
There  are  several  mathematical  difficulties  in  optimizing  for  best  ply 
orientations.  First,  the  objective  function  (total  thickness)  is  not 
directly  a  function  of  angle.  Second,  there  may  be  many  local  minima. 
Third,  if  a  direction  vector  is  found  in  the  combined  angle  and 
thickness  space,  the  magnitude  of  the  scalar  distance  will  have 
different  meaning  for  each  type  of  design  variable.  Finally,  there  is 
the 

the  practical  difficulty  that  ply  orientations  cannot  be  completely 
arbitrary  due  to  manufacturing  limitations.  There  should  be  some 
minimum  angular  step  size  limited  by  the  lay-up  procedures  used.  The 
algorithm  derived  here,  while  not  completly  satisfactory,  attempts  to 
account  for  all  these  difficulties. 

The  approach  taken  is  to  first  divide  the  problem  into  a  multi¬ 
level  optimization  [10],  where  angles  and  ply  ratios  are  optimized 
independently.  We  can  alternate  between  the  two  types  of  optimization 
until  the  laminate  converges  to  a  minimum  thickness  design.  Ply 
thicknesses  are  handled  exactly  as  before,  with  the  given  angles  held 
constant.  During  the  angle  optimization,  ply  ratios  are  held  constant 
and  the  angles  varied  to  minimize  total  thickness. 

The  angle  optimization  used  here  is  not  a  direct  method  like  that 
used  for  the  thicknesses,  but  instead  relies  on  minimizing  a  related, 
unconstrained  function  with  the  assumption  that  total  thickness  will 
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decrease  at  the  sane  time.  One  approach  is  to  chose  a  function  which 
will  lead  to  the  simultaneous  failure  condition,  which  should  result  in 
an  efficient  laminate.  Another  desirable  feature  is  that  the  results 
should  not  be  too  sensitive  to  the  selection  of  initial  angles.  After 
experimenting  with  several  possible  functions,  the  best  was  found  to  be 
the  variance  of  the  all  the  constraints,  given  by  the  equation 
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(27) 


■  © 


*  I 


where  nc  =  m  *  rV 

If  this  function  where  minimized  to  a  value  of  zero,  a  simultaneous 
failure  condition  for  the  laminate  would  be  reached.  In  cases  with 
multiple  loads,  simultaneous  failure  for  all  loads  is  usually 
impossible,  but  we  assume  that  as  the  variance  is  minimized,  as  many 
constraints  will  become  active  as  possible.  It  will  be  shown  that 
simultaneous  failure  is  not  always  the  optimal  condition  for  a  composite 
laminate,  but  for  most  cases  it  will  either  be  the  minimum  or  at  least 
very  close  to  the  minimum  thickness.  Before  trying  the  variance,  the 
author  had  attempted  to  minimize  the  value  of  the  largest  current 
constraint  function.  This  necessitated  finding  some  way  to  handle 
multiple  constraints  that  had  nearly  the  same  value.  This  version  of 
the  program  often  terminated  early  because  a  satisfactory  way  was  never 
derived  for  finding  a  common  vector  that  would  reduce  the  value  of  more 
than  one  constraint  simultaneously.  To  find  the  minimum  of  the 
variance,  a  steepest  descent  method  was  used,  formally,  steepest 
descent  is  considered  the  least  efficient  way  to  minimize  an 
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unconstrained  function,  but  it  was  found  to  be  sufficient  for  the 
current  research.  The  program  should  be  modified  in  the  future  to 
include  a  conjugate  gradient  method  [10]. 

The  steepest  descent,  along  with  most  other  search  methods,  needs 
the  value  of  the  gradient.  Terms  of  the  gradient  are  given  by 
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where 
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i  =  Ze(L)  iG  P  ,£(L)l9i  +  e(L) 1g(P>I ’ei^(L) 
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and 


£(L) >0j  =  ’l  |Al»0j^(L) 

|AU  = 


(29) 


It  should  be  noted  that 

iG(P)|  >0  =  0  for  i  f  P 

The  angular  derivatives  of  the  Q's  and  G's  are  given  in  Appendix  A.  The 
negative  of  the  gradient  will  form  the  search  direction.  The  scalar 
distance  along  the  search  direction  is  found  by  taking  discrete  steps 
and  stopping  when  the  thickness  begins  to  increase  (and  then  taking  one 
step  backwards).  Because  the  variance  is  only  a  function  related  to  the 
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actual  minimum,  we  do  not  determine  the  distance  by  the  magnitude  of  the 
variance,  but  instead,  the  function  we  are  actually  interested  in. 
Thickness  is  calculated  by  using  the  scaling  equations  developed 
previously.  More  efficient  one-dimensional  search  methods  will  have 
difficluties  t«?ith  the  multiple  local  minima. 

The  steps  are  taken  so  that  all  the  angles  change  by  some  minimum 
step.  To  maintain  the  minimum  s.  ep  size,  the  angles  are  incremented  by 
the  equation 

=  0-  +  (CINT  [(k+liZf]  -  CINT  (kZi))A9 

where  CINT  implies  taking  the  closest  integer  value  and  k  is  an 

—V 

incremental  step  counter.  The  direction  vector  Z  is  normalized  by  its 
largest  element.  At  each  unit  increment  of  k,  the  angle  corresponding 
to  the  largest  element  of  3  is  incremented  by  A0.  Other  angles  may  not 
be  incremented  at  each  step,  depending  on  the  relative  values  of  the  Z 
vector  elements.  Thus,  the  direction  vector  is  not  followed  exactly, 
but  rather  on  a  broken  path.  The  amount  of  divergence  from  the  search 
direction  is  determined  by  the  value  of  A9»  If  the  angle  start  out  as 
multiples  of  10°  and  A0  is  10°,  then  the  angles  will  stay  as  multiples 
of  10°  throughout  the  search. 

The  overall  procedure  for  the  multi-level  optimization  can  be 
summerized  as  follows: 

1)  Enter  loads  and  starting  angles 

2)  Find  a  search  direction  based  on  the  variance 

3)  Perform  a  one-dimensional  search  to 
minimize  total  thickness  with  constant  ply 


ratios 


4)  Repeat  from  step  2  until  no  further  changes 
in  angle  can  be  made 

5)  Optimize  the  ply  ratios 

6)  Repeat  from  step  2  until  neither  type  of 
optimization  can  make  further  progress 

Testing  of  the  program  shows  that  one  pass  through  steps  1  to  6  is 
all  that  is  needed.  Usually,  the  angle  optimization  brings  enough  of 
the  constraints  to  critical  values  that  the  ply  ratio  optimization  can 
make  little  progress.  In  turn,  after  the  ply  thickness  routine  is 
finished,  there  is  little  the  angle  optimizer  can  change. 

Typically,  the  angle  optimization  will  need  4-6  search  directions 
to  converge,  requiring  10-20  minutes  for  4  ply  groups  and  a  pair  of 
independent  loads. 


Orthotropic  Laminate 


A  designer  may  not  want  a  general  symmetric  laminate.  He  may  be 
more  comfortable  with  an  orthotropic  laminate  which  eliminates  the  shear 
coupling  terms  and  allows  the  use  of  many  existing  orthotropic  plate 
analysis  equations.  An  orthotropic  laminate  can  be  made  by  keeping  the 
’  aminate  balanced.  That  is,  for  every  ply  at  +8,  there  is  one  at  -8. 
There  may  also  be  manufacturing  reasons  for  wanting  a  balanced  laminate, 
such  as  filament  winding  operations.  There  is  no  difficulty  in 
constraining  the  optimization  procedure  to  yield  ibalanced  laminates. 

Most  sophisticated  optimization  programs  allow  design  variables  to  be 
coupled  so  that  they  maintain  the  same  value.  A  simpler  approach  is  to 
enter  only  the  positive  angle  and  set  the  A^  and  A23  terms  to  zero. 

The  resulting  thickness  found  for  the  positive  angle  must  then  be  split 
between  the  positive  and  negative  angles  in  the  actual  laminate.  With 
the  reduced  A  matrix,  a  faster  matrix  inversion  can  be  written. 

When  designing  with  orthotropic  laminates,  the  orthotropic  axis 
should  not  be  selected  arbitrarily.  For  a  single  load,  the  orthotropic 
axes  should  be  aligned  with  the  principle  axes  of  the  load.  With 
mulitple  loads,  the  selection  is  not  so  obvious.  Finding  the  best  axes 
with  respect  to  the  load  is  a  much  simpler  problem  than  the  general 
optimal  angle  search  discussed  above.  A  search  for  the  best  axes  can  be 
reduced  to  a  one-dimensional  search.  The  procedure  can  be  thought  of  as 
finding  the  best  rigid  body  rotation  of  the  laminate  with  respect  to  the 
loads  while  performing  a  thickness  optimization  of  each  rotation  angle 
(Figure  15).  For  computational  simplicity,  the  program  actually  rotates 


FIGURE  15:  Definition  of  Angles  for  Orthotropic  Laminate 


the  loads  and  keeps  the  laminate  angles  fixed.  Even  this  one¬ 


dimensional  search  could  be  time  consuming  without  a  fast  ply  ratio 
optimization  algorithm.  The  orthotropic  optimization  with  the  strain- 
sphere  failure  criteria  is  fast  enough  to  make  a  search  for  best 
orientation  practical. 

The  search  procedure  can  be  summerized  as  follows: 

1)  Enter  initial  laminate  angles,  loads,  bounds  on 
search  angle,  and  maximum  error  for  search. 

2)  Divide  the  bounded  region  with  4  equally  spaced 
points,  with  endpoints  on  the  bounds 

3)  Find  the  minimum  laminate  thickness  at  each 
point  by  rotating  the  loads  by  the; negative  of 
the  current  test  angle 

4)  Check  for  the  smallest  value  of  the  4 
thicknesses.  The  2  points  on  either  side 
of  the  smallest  one  become  the  new  bounds. 

5)  If  the  bounds  are  greater  than  the  maximum 
error,  repeat  from  step  2.  Only  2  new  points 
need  to  be  calculated. 

The  method  being  used  here  is  very  similar  to  the  bisection  method 
for  finding  the  zero  of  a  function.  Bisection  requires  3  function 
values  in  order  to  reduce  the  size  of  the  region  the  zero  can  be  in. 
Here,  a  fourth  point  is  needed  because  we  are  searching  for  the  zero  of 
the  first  derivative  instead  of  a  zero  of  the  function. 
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III.  APPLICATION 


Examples 


A  few  illustrative  examples  will  be  discussed  to  demonstrate  the 
operation  of  the  optimization  procedures.  A  detailed  comparison  of  the 
weight  savings  possible  with  ply  ratio  optimization,  angle  optimization, 
and  no  optimization  will  be  given  in  the  next  section. 

The  strength  ratios  defined  in  [7]  will  be:  needed  to  show  which 
plies  are  critical  for  given  loads.  The  ratio  is  defined  as  the  value 
of  R  in  the  equation 


R2  e T | 6 [ e  +  RGTe  =  1 

An  R  of  1  means  the  ply  is  at  the  boundary  of  the  failure  envelope.  R's 
1  mean  the  ply  is  in  the  safe  region  on  the  failure  envelope.  The  R's 
can  be  interpreted  as  the  ratio  of  the  applied  load  vector  length  to  the 
maximum  load  vector  length. 

Most  of  the  examples  presented  here  will  use  T300/5208  as  the 
material.  Properties  of  this  material  along  with  Kevlar  and  aluminum 
(used  in  certain  examples)  are  given  in  Table  2.  Figure  16  is  an 
example  output  from  an  Epson  HX-20 

microcomputer.  Only  ply  ratios  are  being  changed  and  the  angles  are 
given  as  0/90/45/-45.  This  example  demonstrates  a  case  where  there  is 
no  severest  load  condition.  Looking  at  the  strength  ratios,  we  can  see 


that  the  90  and  -45  plies  are  near  failure  for  the  first  load  condition, 


Material  Properties 

T300/5298 

E.\=  181  GPa 

EV=  10.3  GPa 

ES=  ?.  1?  GPa 

UX=  .28 

X=  1509  MPa 

X’=  1500  MPa 

V=  40  MPa 

V ’ =  246  MPa 

S=  68  MPa 

Ply  Thickness  .000125  m 

LOADING  1 
N  1=3  MN/m 
M  2-  1  MN/m 
N  6=  O  MN/m 
LOADING  2 
N  1=  1.5  MN/m 
N  2=  1.5  MN/m 
N  6=-.  5  MN/m 


Total  thickness- 
.  0735E-01  in. 

58.  76  Plies 


ANGLE 

RATIO  #PLIES 

0 

.4416  25.95 

90 

.  1236  7.26 

45 

.1774  10.42 

-45 

.2574  15.12 

STRENGTH  RATIOS 

1  ULTIMATE  STRAIN 

>1  IS  S 

AFE 

LOADING 

1 

PLV 

0 

1.4078 

90 

1 

45 

1.  2091 

-45 

1.0973 

LOADING 

2 

PLV 

0 

1 . 0355 

90 

1.  4004 

45 

1. 0331 

-45 

1. 4071 

LAMINATE  STRAINS 
LOADING  1 
el=+3.  628E-03 
e2=+l.  274E-03 
e6*+0.  691E-03 
LOADING  2 
el=+l.  224E-03 
e2=+3.  334E-03 
e6=-2.  157E-03 


Norm.  |  A  | 

in  GPa. 

il06„  20?  I 

20. 0351 

-3.  429 

i  20,  035  i 

51.  6731 

- 1 

-3.  429 

j  -3.4291 

-3. 4291 

— - L_ 

- 1 

24.  309 

Compliance  (normalized) 
in  1/TPa. 


1 

10.  1781 
- - 

-3. 8871 

0.  887 

1 

U- 

-3,  83?  I 
-  | 

21.0201 

2.417 

t 

0.  8871 

2.  41? i 

- - - L. 

- - 1 

41.604 

1 

ENGINEERING  CONSTANTS 

El*  98.  3  GPa 
E2=  47.6  GPa 
E6=  24.  0  GPa 

v21=  0.  382 
v61=  0.087 
vl6=  0.021 


FIGURE  16:  Printout  for  Example  Problem 


and  the  0  and  +45  plies  are  near  failure  for  the  second  load  condition. 
The  normalized  A  matrix  shown  in  the  output  is  defined  as  I aI /h  and  the 
normalized  compliance  matrix  is  the  inverse  of  | A|  times  h. 

The  example  given  in  Figure  17  is  a  case  where  simultaneous  failure 
is  impossible.  The  constraint  curves  in  design  space  are  plotted  to 
show  that  one  constraint  is  never  on  the  boundary  between  the  feasible 
and  infeasible  regions.  The  impossibility  of  simultaneous  failure  is 
also  evident  by  examining  the  failure  envelopes  in  strain  space.  The 
failure  envelopes  for  graphite  epoxy  only  intersect  in  the  first  and 
fourth  quadrants  (Figure  3).  Pure  shear  transformed  to  principle 
strains  is  in  the  second  or  fourth  quardrant.  Even  though  one  ply  is 
never  near  failure,  removing  that  ply  increases  the  total  thickness 
required. 

Table  3  compares  the  results  of  optimization  based  on  the  strain- 
sphere  approximation  and  the  usual  quadratic  interaction  criteria.  The 
ply  ratios  are  quite  close,  demonstrating  that  for  loads  in  the  first 
quadrant,  the  results  are  not  sensitive  to  the  particular  criterion. 
Although  the  approximate  criterion  works  well,  all  results  presented 
elsewhere  in  the  thesis  will  be  based  on  the  quadratic  criteria  unless 
otherwise  stated.  No  detailed  description  of  the  algorithm  for 
optimization  with  the  strain  sphere  is  given,  but  the  method  is  almost 
identical  to  that  used  for  the  quadratic  criteria.  The  major 
differences  are  that  the  gradient  is  redefined  and  the  criterion  only 
needs  to  be  evaluated  for  the  laminate  as  a  whole,  rather  than  for  each 
ply  individually. 

Table  4  is  an  example  of  the  orthotropic  laminate  optimization 
with  optimal  rigid  body  rotation.  The  best  orthotropic  axes  could  not 
have  been  selected  from  inspection  of  the  load  principle  axes.  The 
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n) 


LOAD 
Nl=  0 
N2=  0 

N6=  2  MN/m 


;le 

Ply  Ratio 

#  Plies  Needed 

>5 

.656 

26.9 

>5 

.344 

14.1 

Angle 

45 

-45 


Strength  Ratios 


R 

2.42 

1.00 


Figure  17: 


Constraint  Curves  and  Optimization 
Results  for  45  Under  Pure  Shear 


Nl=  4  MN/m 
N2=  1  MN/m 
N6=  0  MN/ra 


LOADS 


Nl'  =  2.76  MN/m 
N2'=  2.24  MN/m 
N6 ' =-l . 48  MN/m 


Ply  Group 

#  Plies  Needed 
Tsai-Wu 

Approximate 

0 

35.2 

35.2 

90 

7.5 

7.4 

45 

9.9 

10.8 

-45 

33.8 

33.7 

Total  86.5 

87.1 

TABLE  3:  Comparison  of  Approximate  Strain-Sphere  to  Tsai-Wu 

Criteria  for  Optimization 


Two  independent  Loads 


LOADS 

Nl=  2  MN/m 

Nl'=  1.25  MN/m 

N2=  1  MN/m 

N2'=  1.75  MN/m 

N6=  0  MN/m 

N6 1 =  -.43  MN/m 

# 

Plies  Needed 

Angle 

Fixed  Axis 

Variable  Ortho.  Axis 

0 

17.2 

11.2 

90 

17.2 

3.9 

45 

8.6 

16.3 

+45 

8.6 

16.3 

Total 

51.6 

47.7 

TABLE  4:  Comparison  of  Optimal  Orthotropic  Laminates 

with  Fixed  and  Variable  Orthotropic  Axis. 

Optimal  orthotropic  axis  at  -3<?  . 


results  for  an  orthotropic  laminate  with  the  axes  arbitrarily  set  on  one 
of  the  load  principle  axes  are  also  given.  The  difference  is 
substantial.  Both  examples  are  based  on  the  maximum  strain-sphere 
criterion. 

Angle  optimization  is  only  needed  if  there  is  more  than  one 
independent  load.  For  a  single  load,  the  algorithm  will  simply  rotate 
the  plies  so  that  they  lie  on  the  load  principle  axes.  This 
characteristic  shows  that  there  is  more  than  one  minima,  since  an  angle- 
ply  (consisting  of  a  48  and  a  -0  ply  group)  is  more  efficient  than  a 
cross-ply  laminate  (0's  and  90' s  only).  The  program  does  not  converge 
to  the  angle-ply  solution  unless  the  initial  angles  are  close  to  the 
final  value.  We  cannot  predict  the  result  when  multiple  loads  are 
included.  To  show  the  relationship  between  load  principle  axis  and 
optimized  ply  orientations,  2  independent  loads  that  fall  on  the  same 
Mohr's  circle  have  been  used  as  the  design  requirements.  The  loads  and 
ply  orientations  can  be  superimposed  on  the  same  Mohr's  circle.  Figure 
18  reflects  some  of  the  symmetries  of  the  optimized  laminate.  An 
interesting  example  of  how  non-intuitive  composites  can  be  is  shown  in 
Figure  19.  Two  equal  magnitute  uniaxial  loads  <.  re  entered  with  one  of 
the  loads  rotated  by  -40  from  the  laminate  axis.  Instead  Ox  placing  the 
plies  on  the  principle  axes,  the  computer  has  picked  slightly  different 
angles,  which  give  a  thinner  laminate  than  if  the  principle  axes  had 
been  used.  The  starting  angles  were  0/90/45/-45,  but  the  angles  have 
converged  so  that  only  2  ply  groups  remain. 

Although  there  is  now  a  method  for  finding  good  ply  orientations, 
we  still  need  to  know  how  many  initial  angles  should  be  used,  and  their 
initial  values.  One  reason  the  search  based  on  constraint  variance  was 
selected  is  because  it  seems  to  be  less  sensitive  to  choice  of  initial 


PLY  FIBER  DIRECTIONS 


FIGURE  18 :  Mohr's  Circle  Representation  of  2  Independent  Loads 
with  Superimposed  Optimized  Ply  Orientations 


Initial  angles  (0/90/^45) 
All  angles  plotted  as  20 


y  ply  fiber  direction 


FIGURE  19:  Mohr's  Circle  Representation  of  2  Independent  Load 
with  Superimposed  Optimized  Ply  Orientations 


Initial  angles  (0/90/+45) 
All  angles  plotted  as  20 


angles  than  some  of  the  other  methods  tried.  The  number  of  angles 
needed  is  still  an  open  question.  A  quick  look  at  gradient  information 
suggests  that  too  few  angles  (2  for  example)  will  make  the  laminate 
sensitive  to  small  changes  in  orientation  or  load.  The  0/90/45/-45 
starting  point  selected  for  all  the  above  examples  has  been  found  to 
give  effiecient  laminates  without  the  complexity  of  adding  a  lot  of 
angles.  Most  of  the  examples  run  where  with  2  loads,  but  a  couple  of 
cases  were  tried  with  4  loads.  The  4  ply  group  laminate  was  still 
adequate  despite  the  additional  loads. 

All  the  examples  given  were  run  by  applying  the  angle  optimization 
first  and  then  the  ply  ratio  optimization.  After  the  ply  ratio 
optimization,  no  further  attempt  at  changing  the  angles  was  made. 

There  is  the  possiblility  that  the  combined  angle/ply  ratio 
optimization  will  yield  a  laminate  with  total  thickness  greater  than 
would  have  been  produced  by  ply  ratio  optimization  alone.  By  bringing 
more  constraints  into  play,  the  angle  optimization  may  prevent  the  ply 
ratio  program  from  making  as  much  progress  as  it  would  have  starting 
from  some  arbitrary  initial  angles.  Often,  the  ply  ratio  program  will 
not  be  able  to  change  the  laminate  at  all,  leaving  the  ply  ratios  equal. 
From  the  evaluation  presented  later  in  this  thesis,  we  can  see  that 
there  is  a  choice  of  which  variables  are  optimized.  There  may  be  some 
motivation  for  keeping  the  ply  ratios  constant,  or  near  constant.  In 
which  case,  angle  optimization  will  still  give  an  efficient  laminate. 

If  angles  are  fixed,  ply  ratio  optimization  alone  will  also  give  an 
efficient  laminate. 

The  capability  to  optimize  hybrid  laminates  is  easily  added  to  the 
existing  programs.  When  the  A  matrix  is  formed,  the  Q’s  associated  with 


the  proper  material  are  used.  Also,  the  constraint  test  and  gradient 
calculations  must  use  the  appropriate  values  of  the  G's  for  whichever 
material  the  given  ply  is  made  from.  The  example  given  in  Table  5 
shows  the  results  for  a  hybrid  made  from  alternating  ply  groups  of 
graphite/epoxy  and  Kevlar/epoxy,  with  each  orientation  duplicated  by 
both  materials.  For  strength  constraints,  the  Kevlar  is  usually 
completly  removed.  The  combination  of  glass/epoxy  and  graphite/epoxy 
was  found  to  give  similar  results.  No  strength  advantage  has  been  found 
by  going  to  hybrid  systems. 


Nl=  4  MN/m 
N2=  1  MN/m 
N6=  0  MN/m 


Material 

Angle 

#  Plies  Needed 

Graphite 

0 

42.0 

Kevlar 

0 

0.0 

Graphite 

90 

4.9 

Kevlar 

90 

0.0 

Graphite 

45 

9.4 

Kevlar 

45 

0.0 

Graphite 

-45 

9.4 

Kevlar 

-45 

0.0 

TABLE  5:  Hybrid  Laminate  Example 


Potential  Weight  Savings 


Optimization  would  be  of  little  interest  if  the  potential  gains 
were  only  a  few  percent.  In  fact,  for  strength  controlled  laminates, 
the  weight  savings  are  usually  in  the  range  of  20-50%,  as  compared  to 
quasi-isotropic  lay-ups.  The  thickness  difference  due  to  optimization 
with  a  single  biaxial  load  can  be  seen  in  Figure  20.  This  is  a  fairly 
general  graph,  since  any  biaxial  load  can  be  transformed  to  a  shear-free 
axis  (principle  directions)  and  differences  in  N1  would  just  cause  a 
proportional  change  in  total  thickness.  It’s  interesting  to  note  that 
the  0/90/45/-45  laminate  is  thinner  than  the  0/90.  Beyond  a  load  ratio 
of  about  2  (N  ]/N  2),  the  9CP  ply  in  the  0/90/+45/-45  laminate  is  dropped 
completly,  making  a  tri-directional  laminate  that  is  more  efficient  than 
the  0/90.  A  good  rule  in  design  is  to  make  the  laminate  axes  and  load 
principle  axes  coincide  when  there  is  a  only  a  single  load.  The  angle 
optimization  routine  will  give  this  intuitive  result.  However,  with  4 
or  more  available  orientations,  the  ply  ratio  optimization  is  forgiving 
if  the  principle  directions  are  not  used.  A  0/90/45/-45  laminate  was 
rotated  as  a  rigid  body  with  respect  to  a  fixed  4:1  biaxial  load.  The 
laminate  was  optimized  at  5°  increments  of  rotation.  The  difference 
beween  the  thickest  and  thinnest  resulting  laminate  was  only  5%. 

When  two  or  more  independent  loads  are  combined,  the  anisotropic 
advantage  of  composites  becomes  less  significant,  (because  there  is  less 
of  a  distinct  preferred  direction)  but  the  savings  due  to  optimization 
can  still  be  substantial.  Because  tbnre  are  an  infinite  number  of  load 
combinations,  it's  impossible  to  draw  any  general  graphs  demonstrating 
the  gains  due  to  optimization.  To  give  an  indication  of  the  trends,  a 
series  of  18  load  combinations  was  devised,  where  each  load  combination 


consists  of  a  pair  of  biaxial  loads.  Because  of  the  directionality  of 
composites,  loads  with  differing  principle  axes  are  of  greatest  interest 
for  excerising  the  procedure.  The  load  combinations  and  principle  axes 
orientations  are  given  in  Figure  21.  This  group  of  load  cases  is  not 
intended  to  be  all-encompassing,  but  represents  some  worst  case 
conditions  for  taking  advantage  of  a  directional  material.  Most  of  the 
loads  are  in  tension,  although  cases  13-15  are  compression-compression 
and  cases  16-18  are  mixed  tension  and  compression.  The  magnitudes  of 
the  principle  components  of  the  loads  have  been  made  equal  in  most  of 
the  cases  in  order  to  ensure  both  loads  influence  the  final  design.  A 
small  load  might  never  form  part  of  the  boundary  between  feasible  and 
infeasible  design  space.  Intial  angles  are  0/90/45/-45  for  all  the 
types  of  optimization  considered  below.  The  next  section  will  show  that 
equal  angular  spacing  is  a  good  starting  point  for  picking  angles  for 
the  optimization  code  to  work  with.  Ply  ratio  optimization 

alone  will  be  considered  first.  Figure  22  shows  the  weight  savings  of 
optimized  0/90/45/-45  laminates  versus  unoptimized  laminates  of  the  same 
angles.  Kevlar  material  was  taken.  Again,  the  load  cases  are 
arbitrary,  but  the  point  to  be  made  is  that  around  a  25%  weight  savings 
can  be  expected  from  using  optimization  with  a  wide  variaty  of  loads. 

In  some  cases  the  savings  can  be  even  larger  (40-50%  for  several  of  the 
load  cases).  To  show  that  the  results  are  not  material  dependent,  the 
same  loads  have  been  applied  to  laminates  made  from  graphite/epoxy 
(T300/5208).  This  time  the  savings  are  compared  to  aluninun,  (with 
density  diffei  nee  included).  The  large  differences  between  the 
optimized  and  unoptimized  laminates  are  still  evident  (Figure  23). 

The  first  12  load  cases  (all  tension-tension  loads)  were  used  to 
test  the  strain-sphere  criterion.  When  averaged  over  the  12  loads,  this 
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FIGURE  21:  Load  Case-  Matrix  for 


Independent  Loads 
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FIGURE  22:  Weight  Savings  of  Optimized  Kevlar  Laminates  Over  Quasi-Isotropic 

Load  cases  correspond  to  Figure  21 


EQUAL  PLY  RATIO 


O 


approximate  criterion  was  found  to  be  only  7 Z  conservative  as  compared 
to  the  quadratic  criterion.  Thus,  when  only  tension  loads  are 
considered  (or  with  small  compression  components),  the  approximation  nay 
be  desirable  if  computation  time  is  a  factor. 

The  orthotropic  axis  optimization  is  based  on  the  strain-sphere 
criterion.  This  type  of  optimization  was  also  tested  against  the  first 
12  load  cases.  The  results  are  presented  in  Table  6.  The  average 
thickness  is  nearly  the  same  as  for  ply  ratio  optimization  alone, 
despite  the  conservative  criterion  and  the  added  constraint  of 
maintaining  a  balanced  laminate. 

Finally,  angle  optimization  was  also  applied  to  laminates  subjected 
to  all  18  load  cases,  both  with  and  without  subsequent  ply  ratio 
optimization.  A  minimum  angle  change  of  5  was  always  taken  (see 
equation  31).  With  angle  and  ply  ratio  optimization,  the  average  weight 
savings  is  about  6.5%  better  than  ply  ratio  optimization  alone,  but  the 
results  for  individual  cases  vary  widely.  Some  load  cases  resulted  in 
slightly  greater  thickness  with  angle  optimization  than  without.  The 
results  are  almost  identical  if  angle  optimization  is  used  without  ply 
ratio  optimization  at  all.  This  demonstrates  that  the  2  types  of  design 
variables  are  almost  redundant,  and  optimizing  both  is  usually  not 
required. 

As  Table  6  demonstrates,  the  designer  has  some  options  for  picking 
the  parameters  to  be  optimized.  The  final  results  do  not.  vary  much  for 
either  ply  ratio  optimization  with  fixed  orientations,  orthotropic 
laminates  (with  rigid-body  rotation  allowed),  or  angle  optimization 
alone.  The  degree  of  strength  anisotropy  appropriate  to  the  design 
requirements  can  be  achieved  by  varying  any  of  these  parameters.  This 
means  that  composite  materials  have  even  more  flexibility  than 
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previously  imagined.  The  parameters  to  be  optimized  can  be  constrained 
by  other  considerations  (such  as  manufacturing)  and  efficient  laminates 
can  still  be  produced. 
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Number  of  Angles  Necessary 


To  use  any  of  the  methods  described  in  this  thesis ,  the  number  of 

ply  orientations  initially  given  to  the  optimization  program  must  be 

chosen.  The  performance  of  various  laminates  with  different  numbers  of 

initial  angles  was  investigated  to  give  some  indication  of  how  to  pick 

these  angles.  A  likely  starting  point  for  initial  angles  is  to  space 

the  ply  angles  evenly  over  the  180  available.  This  class  of  laminates 

will  be  referred  to  as  tt /n  laminates,  where  n  is  the  number  of 

orientations  in  the  lamiante.  A  tt/4  laminate  has  an  angular  spacing 

o 

between  ply  groups  of  45.  These  lamiantes  are  quasi-isotropic  for  n 
greater  than  2  [7],  This  is  a  reasonable  starting  point  for 
optimization  since  there  are  no  preferred  directions  to  initially  bias 
the  result. 

The  total  thickness  turns  out  to  be  almost  independent  of  the 
number  of  angles  for  a  single  biaxial  load  (Table  7).  By  applying  the 
18  load  cases  introduced  in  the  last  section,  a  comparison  for  multiple 
loads  can  also  be  made.  The  average  weight  savings  (compared  to  a  0i/90i 
/45^/-45^  without  optimization)  is  given  in  Table  8.  For  n  greater  than 
3,  the  averages  are  very  close.  It  is  a  little  deceptive  to  take  the 
average.  When  examined  case-by-case,  the  thickness  differences  between 
the  types  of  laminates  can  be  great  for  a  particular  load  case  (Figure 
24).  These  differnces  may  be  largely  due  to  numerical  problems.  With 
a  large  number  of  ply  groups,  the  program  may  occasionally  terminate 
early  because  of  the  large  number  of  simultaneously  active  constraints. 
Despite  this  variation,  the  7r/4  laminate  seems  to  be  adequate  for 
multiple  loads.  Increasing  the  number  of  angles  will  not  guarantee  a 
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Nl=  3  MN/m 
N2=  1  M/m 
N6=  0  MN/m 


#  Ply  Groups  Total  #  of  Plies 


o 

3 

52 

45 

4 

49 

30 

6 

51 

18 

10 

50 

10 

18 

51 

Table  7:  Total  Thickness  Required  to  Support  a 
Single  Load  for  Various  / n  Laminates 


i^Ply  Groups  %  Weigth  Savings 


60 

3 

19 

45 

4 

23 

30 

6 

24 

10 

18 

25 

TABLE  8:  Average  Percent  Weight  Savings  Over  Quasi- 

Isotropic  for  All  13  Combined  Load  Cases 


FIGURE  24:  Weight  Savings  of  Various  71/n  Laminates 


LOAD  CASE 


better  laminate 


The  examples  in  this  study  included  some  tt/  1 8  laminates.  An  early 
idea  was  to  find  optimal  angles  by  looking  at  a  large  number  of  initial 
angles  and  seeing  what  remained  after  ply  ratio  optimization.  The 
actual  result  is  a  little  surprising.  Instead  of  a  few  optimal  angles 
dominating  the  f'nal  laminate,  the  ply  ratios  plotted  against  angle  form 
almost  a  continous  function  (Figure  25).  All  18  ply  groups  are  near 
failure  for  this  laminate.  For  some  multiple  load  test  cases,  2  peaks 
in  this  pseudo-continous  function  would  form.  A  case  that  formed  more 
than  2  peaks  was  never  found. 

In  conclusion,  the  number  of  initial  angles  can  be  bounded  to  a  few 
choices.  With  only  2  orientations,  we  must  have  some  way  of  picking  the 
angles  since  the  unoptimized  laminate  will  have  a  directional 
preference.  There  doesn't  seem  to  be  any  advantage  to  using  more  than  4 
orientations.  Thus,  n/A  lamiantes  were  used  for  most  of  the  examples  in 
this  thesis,  and  are  suggested  as  a  starting  point  for  design. 


IV.  ANALYTIC  STUDIES 


Maximum  Strain  Energy  Density 


A  visual  representation  of  how  a  laminate  adapts  to  the  given  load 
requirements  would  be  desirable.  A  conventional  failure  envelope 
representation  is  not  acceptable  because  with  multiple  loads,  3- 
dimensions  would  have  to  be  shown  in  order  to  account  for  the 
differences  in  shear  between  the  loads.  The  approach  taken  here  is  to 
plot  the  maximum  strain  energy  density  the  laminate  can  sustain  as  a 
function  of  load  principle  axes  orientation  with  respect  to  the  laminate 
axes.  Then,  on  the  same  graph,  the  strain  energy  density  actually 
produced  by  various  loads  (in  particular,  the  design  loads  )  can  also  be 
plotted.  There  is  a  loss  of  information  in  such  a  graph.  The 
combination  of  N^.  to  N  (magnitutes  of  loads  on  the  principle  axes) 
that  produces  the  maximum  strain  energy  is  an  intermediate  calculation 
and  would  not  be  displayed.  The  graph  is  not  really  a  failure 
representation,  since  it  would  be  possible  to  have  loads  which  produced 
less  strain  energy  but  still  caused  failure.  Despite  these  limitations, 
these  graphs  do  give  a  good  intuitive  feel  for  the  characteristics  of  an 
optimized  laminate. 

The  approximate  strain-sphere  failure  criterion  is  the  starting 
point  for  the  derivation.  We  assume  the  maximum  strain  energy  occurs 
when  the  failure  criterion  reaches  an  equality.  Then 


1  E2 
2  e6 


(33) 


There  are  no  shear  loads,  so  that  the  stress-strain  relation  can  be 


76 


written 


(e)  =  N-,  [A-1  J 


1 

A 

0 


where  A  is  defined  by 


A  =  N2/N1 

The  average,  laminate  strain  energy  density  is  given  by 

U  =  TFr  fe>TU|  (d  i 

where  h  is  the  total  thickness.  Substituting  equation  (34)  into  (36) 

yields 
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2h  (an  +  2a12A  +  a22A  ^ 


(■ 


where  a^'s  are  elements  of  the  inverted  A  matrix. 
Substituting  equation  (34)  into  the  failure  criterion  yields 


[ (a-|  ■  +  a^A)2  +  (ai 2  +  a22A^2  +  \  ^a13  +  a23A^  =  1)2 

N]2  =  b2/C(af1  +  a122  +  \  ai32)  +  (2a11a12 

+  2a19a99  +  a19a99)A  +  (a  2  +  a  2  +  l  a  2)A2]  (38) 


Q  =  2{i  a 1 3a23  +  alla12  +  a 1 2a 22 


R  =  (a122  +  a222  +  \  a232) 


=  b2/(P  +  QX  +  RX2)  (4 

Substituting  (40)  into  (38),  energy  density  becomes 

2 

all  +  2a12x  *  a22X 

~  2h  P  +  QX  +  RX2  (4 

A  derivative  with  respect  to  X  is  taken  in  order  to  find  the  maximum 


value . 


dX  ~  2h  ^2a12  +  2a??^ )(P  +  QX  +  RX2) 


12  T  22 


-  (an  +  2a12X  +  a2  X2)(Q  +  2RX)]/(P  +  Qx  +  RX2)2 


If  du/dX  =0,  then 


(2a12  +  2a22A)(P  +  Q*  +  R^2)  "  (*n  +  2a]2X2)(Q  +  2RX)  =  0 


322Q  '  2a12R 
2( ^22P  ~  a-|  i  R) 

2a12P  '  all^ 


By  substituting  (43)  into  (42),  X  can  be  written  as 

i  _  -B  ±  B2  -  4AC 
Amax  - 2A - 

and  the  maximum  strain  energy  density  is  given  by 


max 


b2 

2F 


ln 


+  2al?Xm^v  +  aooA 

\c  max  22  max 


max  max 


(45) 


To  use  the  relations,  a  rigid  body  rotation  is  performed  on  the 
laminate,  changing  all  the  angles  by  the  rigid  body  rotation  angle  f  . 
The  A  matrix  and  its  inverse  is  calculated  for  the  the  new  angles. 

Both  values  of  A  are  substituted  into  equation  (45)  and  the  one 
yielding  the  largest  strain  energy  is  taken. 

Figure  26  shows  a  typical  graph  for  an  optimized  laminate.  The 
strain  energy  density  actually  produced  by  the  design  loads  are  also 
plotted  as  points.  We  can  see  how  the  laminace  has  adapted  to  these 
loads.  The  function  has  to  repeat  after  90#because  in  the  derivation,  N 
and  Nix  are  interchangeable.  In  Figures  27-29  the  graphs  are  for 
laminates  optimized  to  a  pair  of  loads  with  equal  principle  magnitutes 
but  with  different  angular  spacings  between  their  principle  axes  of  the 
loads.  The  graphs  show  that  as  the  angular  spacing  increases,  the 
laminate’s  degree  of  aniostropy.  decreases.  If  there  are  many  loads  of 
near  equal  magnitude,  and  with  widely  spaced  principle  axes,  then  the 


laminate  would  have  to  be  quasi-isotropic.  There  is  a  limit  to  how 
adaptable  the  laminate  can  be.  The  strain  energy  density  will  be  close 
to  a  sin  48  function,  no  matter  how  many  ply  groups  are  available. 


Design  Principle  Loads  and  Orientations 

N,=  2  MN/m  N ' =  4  MN/m 

A  i 

NIl;=lMM/m  N|I=1  MN/m 

\//  =  0°  \p=  40° 


FIGURE  26:  Strain  Energy  Density  Versus  Principle  Direction 
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Optimality  Criterion 


>1 
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The  question  of  what  constitutes  an  optimized  laminate  (besides  the 
statement  that  it  has  minimum  thickness)  can  be  approached  by 
considering  what  equality  conditions  must  be  true  at  the  optimum.  This 
is  called  an  optimality  criterion  approach.  Some  existing  optimization 
programs  [3]  are  based  on  the  assumption  that  strain  energy  density  will 
be  equal  in  all  the  plies  at  the  optimum.  This  kind  of  criterion  is 
based  on  experience  with  other  types  of  structures,  such  as  trusses. 

The  failure  criterion  doesn't  influence  the  selection  of  ply  ratios,  but 
only  the  total  thickness  scaling. 

The  strain-sphere  criterion  is  simple  enough  that  for  single 
loading  conditions,  an  optimality  criterion  can  be  derived  directly  from 
the  failure  equation.  Taking  only  ply  group  thickness  as  the  design 
variables,  the  minimum  thickness  point  can  be  found  from  the  Langrange 
multiplier  equation 


Vh  +  XVC  =  0 


Terms  of  the  gradient  of  the  constraint  can  be  written  as 


2 1  |  t|  e  h 
1  '  »  hi 


where 


1  0  0 

T  -  0  1  0 
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From  equation  (21),  we  have  the  result  that 
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Substituting  into  the  Langrange  multiplier  equation  (46),  for  each 
component,  we  have 


,-*T 


1  +  |T|  e>h.  =  0 


(50) 


Thus,  each  ply  group  must  satisfy  the  equation 


eT|T|  |A-1 |  |Q(j)|  t-l 


(51) 


where  ^  is  the  same  constant  for  each  ply  group. 

The  strain  energy  density  criterion  could  be  written  as 


cT|Q<1>|  t  -  X 


(52) 


which,  again,  must  be  satisfied  for  each  ply  group.  There  is  a 
significant  difference  between  the  two  criteria.  The  implications  of 
equation  (51)  should  be  studied  in  more  detail.  Perhaps  a  more  direct 
solution  to  the  optimization  problem  can  be  found. 
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V.  CONCLUSIONS 


A  series  of  effective  laminate  optimization  programs  have  been 
developed  and  thoroughly  tested.  The  programs  have  been  designed  to  be 
compact  and  efficient  enough  to  operate  on  the  some  of  the  smallest 
microcomputers.  Although  not  as  general  or  sophisticated  as  some  of  the 
optimization  codes  currently  available,  these  programs  offer  good 
performance  and  are  very  easy  to  use  even  for  those  unversed  in 
optimization.  No  program  in  the  literature  has  been  found  that  can 
perform  angle  optimization  or  the  orthotropic  axis  optimization.  Thus, 
much  greater  flexibilty  is  now  available  to  the  designer. 

The  gains  due  to  optimization  have  been  found  to  be  substantial, 
with  typically  a  30%  weight  savings  as  compared  to  quasi-isotropic 
laminates.  Surprisingly,  these  large  gains  can  be  made  with  either  of  a 
couple  of  design  parameters.  The  designer  can  either  optimize  the  ply 
ratios,  or  the  angles  and  usually  get  equally  efficient  laminates.  Or, 
he  may  chose  to  constrain  the  laminate  to  be  orthotropic  after 
optimization.  If  the  orthotropic  axis  is  free  to  change,  efficient 
laminates  can  be  designed. 

By  trying  many  example  cases,  it  has  been  found  that  a  tt/4  laminate 
is  a  good  starting  laminate.  By  starting  with  quasi-istropic  laminates, 
no  knowledge  of  desired  starting  orientations  for  the  particular  loads 
is  needed.  Increasing  the  number  of  initial  orientations  does  not  seem 
to  improve  the  final  laminates. 

An  approximate  failure  criteria  has  been  found  to  give  good  results 
while  substantially  decreasing  the  computation  times  needed.  The 
approximate  criteria  could  be  particularly  important  when  the 
optimization  procedure  is  tied  into  a  finite  element  code  on  an 
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iterative  basis,  where  the  repeated  optimizations  could  become 
excessively  time  consuming. 

The  approximate  criteria  also  allows  some  analytic  studies  of 
optimized  laminates.  A  representation  of  the  optimized  laminates 
strength  anisotropy  has  been  developed  based  on  the  maximum  strain 
energy  density.  Graphs  made  with  this  formulation  show  how  the 
laminates  match  the  load  requirements.  Also,  there  is  a  limit  to  the 
adaptablility  of  a  laminate.  As  more  load  requirement  are  added, 
eventually  the  laminate  must  become  quasi-isotropic.  An  optimality 
criterion  can  also  be  derived  from  the  approximate  failure  criterion 
which  can  be  the  subject  of  future  investigations. 

Hopefully,  tailored  laminates  will  come  more  common  as  these  new 
tools  are  made  available  to  designers,  enhancing  the  desirability  of 


composites 
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APPENDIX  A 


v. 


Angular  Derivatives 

The  derivatives  of  the  stiffness  and  failure  parameters  are  found  using 
the  multiple  angle  transformations  of  Tsai  [5] 
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3Q 

30 
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30 


—  =  -2U2  sin  20  -  4U3  sin  40 

—  =  2U2  sin  20  -  4U3  sin  40 

—  =  4U3  sin  40 

—  =  4U3  sin  40 

—  *  U2  cos  20  +  4U3  cos  40 

—  =  U2  cos  20  -  4U3  cos  40 


where 

u2  -  «xx  -  V 

u3  '  «xx  +  V  '  2<>xy  '  «ss> 

Partials  of  6^  can  be  found  with  the  same  equations,  but  with 

U2  *  >'2  <Gxx  -  V 

U3  *  W  <Gxx  +  Gyy  *  2Gxy  '  4Gss> 

The  linear  terms  of  the  failure  equation  has  become 

3G1 

3Q-  *  -2q  sin  20 
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3Gn 

gg-  =  2q  sin  29 


■  =  2Q  cos  29 


where 


APPENDIX  B 


Program  for  Ply  Ratio  and  Angle  Optimization 

The  following  program  optimizes  composite  laminates  for  minimum 
weight  subject  to  inplane  strength  requirements.  Program  options  are: 

1)  optimized  ply  ratios  2)  optimize  ply  angles  and  ratios  3)  perform 
laminate  plate  analysis  without  optimization.  Inputs  include  initial 
ply  angles,  loads  (multiple  independent  loads  possible)  and  a  material 
selection.  Material  properties  for  common  composites  are  stored  in  a 
library,  or  new  properties  can  be  entered  by  following  prompts.  The 
program  is  interactive  and  use  should  be  obvious  from  displayed  prompts. 
A  typical  computer/user  dialogue  is  given  below,  along  with  the 
resulting  output. 

The  program  is  written  for  an  Epson  HX-20  microcomputer  which  uses 
a  fairly  standard  form  of  BASIC.  The  major  exception  are  the  GET%  and 
PUT%  commands  to  addrss  the  material  library.  These  can  be  replaced  by 
disk  file  operations  on  most  other  computers.  The  other  possible  change 
would  be  the  explicit  double  precision  symbol  used  in  the  program. 
Test  have  shown  that  double  precision  is  not  really  needed  and  could  be 
left  out  when  using  other  machines. 


COMPUTER/USER  DIALOGUE 


LCD  Display 


Keyboard  Response 
(comments  in  parenthesis) 


Press  any  key  when  desired 
Material  appears 

T300/5208 

B  (4) /5505 

AS/3501 

Scotchply  1002 

Kevlar  49/Epoxy 

Aluminum 

New 

REVIEW  OR  NEW  DATA  (R/N)  ? 

WHICH  MATERIAL  WILL  YOU 
REPLACE  (0-5)  ? 

EX(GPa)  =  ? 

EY(GPa)  =  ? 

VX  =  ? 

ES(GPa)  =  ? 

X (MPa)  =  ? 

X' (MPa)  =  ? 

Y (MPa)  =  ? 

Y' (MPa)  =  ? 

S (MPa)  -  ? 

THICKNESS  (m. )  =  ? 

NAME  (15  CHR  MAX)  ? 

ADDITIONAL  CHANGES  ? 


RUN  RETURN  (unless  otherwise 
noted,  "Return"  key  pressed 
after  each  keyboard  entry) 


(random  key  pressed  when 
"New"  appears  on  screen) 

N 

5 

(materials  numbered  in  same 
order  as  listed  T300/5208=0) 

185 
6.76 

.2 

5.86 

680 

690 

(primed  constants  imply 
compressive  properties) 

16 

186 
72 

125E-6  (ply  thickness) 
HMS/3002M 

N 


LCD  Display 


Keyboard  Response 
(comments  in  parenthesis) 


HOW  MANY  PLY  GROUPS  ? 
ENTER  PLY  GROUP 
ORIENTATIONS 
PLY  1  =  ? 

PLY  2  =  ? 

PLY  3  =  ? 

PLY  4  =  ? 

ENTER  NUMBER  OF  INDEPEN¬ 
DENT  LOADING  CONDITIONS  ? 

LOAD  1  in  MPa 

Nl  =  ? 

N2  =  ? 

N6  =  ? 

LOAD  2  in  MPa 
Nl  =  ? 

N2  =  ? 

N6  =  ? 

OPTIMIZATION  OR 
ANALYSIS  (0/A) 

RATION  OR  ANGLE 
OPTIMIZATION  (R/A) 


(Materials  list  begins  again, 
this  time  with  the  new  material 
replacing  aluminum,  when  it 
appears  a  key  is  pressed) 

4 

0 

90 

45 

-45 

2 

3 

2 

.5 

1 

4 

0 

0 

R 


WORKING  ITERATION  1 

TOTAL  THICKNESS  = 

1.71342  E-02  m. 

137.07  PLIES 
HIT  ANY  KEY  TO  CONTINUE 

Press  Y  if  printout  of  displayed 
result  is  desired.  Press  N  if 
not 

PLY  PROPERTIES 


LOADS 

TOTAL  THICKNESS  & 
PLY  RATIOS 

STRENGTH 

LAMINATE  STRAINS 
STIFFNESS  MATRIX 
COMPLIANCE  MATRIX 
PLY  RATIO  GRAPH 


FINISHED 

HIT  ANY  KEY  TO 

CONTINUE 


(after  4  iterations  and  about 
7  minutes  the  computer  beeps 
that  the  solution  has  been 
found.  This  example  ran  for 
an  unusually  long  time.  Most 
problems  will  run  in  less  time) 

(press  any  key,  no  return) 


Y  (return  key  not 
used  for  these  responses) 

Y 


Y 

Y 

Y 

Y 

Y 

Y 

(after  entire  list  of  print 
out  options  is  presented, 
computer  produces  the  print 
out  shown  on  next  page) 

(pressing  a  key  restarts 
program.  Press  "BREAK" 
key  to  exit) . 


wmmmmmmmmmmm 

Material  Properties 

HMS/3002M 

EX=  185  GPa 

EV=  6.  76  GP3 

ES=  5.  86  GPa 

UX=  .2 

X®  680  MPa 

X’  =  690  MPa 

V=  16  MPa 

V'  =  186  MPa 

8=  72  MPa 

Ply  Thickness  .000125  m 

LOADING  1 
N  1=  3  MN/m 
N  2=  2  MN/m 
N  6=  .  5  MN/m 
LOADING  2 
N  1=  1  MN/m 
N  2=  4  MN/m 
N  6=  0  MN/m 


Total  thickness® 
.  0171E+00  m. 

137.  07  Plies 


ANGLE 

RATIO 

#PLIE: 

0 

.  3476 

47.  65 

90 

.  5281 

72.  38 

45 

.  1243 

17.  04 

-45 

0 

0 

STRENGTH  RATIOS 
1=ULTIMATE  STRAIN 
>1  IS  SAFE 
LOADING  1 


PLV 

0 

1.  1622 

90 

1 

45 

1.  1616 

LOADING 

2 

PLV 

0 

1.0053 

90 

1.  2115 

45 

1.  0122 

LAMINATE  STRAINS 
LOADING  1 
el=+2.  182E-03 
e2=+0.  902E-03 
e6=+l.  096E-03 
LOADING  2 
el=+0.  697E-03 
e2=+2.  216E-03 
e6=-l.  467E-03 


Norm.  |A| 

in  GPa. 

i  74.7621 

6.  518 i 

5.  548 

'i  6.5101106.9671 

5.  548 

i  5. 5481 
>  1 

5.  548 i 

11.  016 

Compliance 

''normalized) 

in  1/TPa. 

i  13.  921 ! 

-0.  497 i 

-6. 760 

1  1 

I  -0.  497  i 

9.6171 

-4. 593 

-| 

i  -6.  760  i 
‘  1 _ 

-4.  593 j 

_ _ _ L_ 

96.  497 

_ i 

ENGINEERING  CONSTANTS 

El=  71.8  GPa 
E2=184.  0  GPa 
E6=  10.4  GPa 
v21®  0.036 
v61=-0.  486 
ol6=-0.  070 


Output  Produced  from  Example  Dialogue 


Reproduced  Irom 
beil  available  copy. 


10  '**  MAIN  CLASS** 

20  CLEAR  75,330 

30  WIDTH  20,4 

40  DEFFIL  55,0 

50  DEFINT  I-PSDEFDBL  F 

60  DIM  A<3»3),B<6>9),C<6 

,6>,D<3,3>,G<3,3),XN<4,3 

>,AI<3,3),Q<3»3),H<6),R< 

3>,S<3>,T<6>,U<5),U<7>,X 

<6>,V<3),Z<6),E(4,3) 

65  DIM  W(24>  6) , CON <24 ) 

70  DIM  Ck<10,2> 

80  DEF  FNDEG<X)=X*57.  295 
78 

90  DEF  FNRAD<X)=X/57.  295 
78 

100  ’**  MAIN  ** 

105  RESTORE 

110  READ  IMAX,E2,E5,E6 

120  ITER=1 

130  GOSUB  2540 

140  CLS: PRINT  "OPTIMIZAT 

ION  OR": INPUT  "ANALYSIS 

<0/A)"5A$ 

150  IF  A$="A"  THEN  6500 
152  INPUT"RATIO  OR  ANGLE 
OPT  <R/A)"!A$:IF  A*="A 
"  THEN  INPUT  "DELTA"; DEL 
TA 

155  CLS:  PR I NT"WORK I NG " : 
PRINT" ITERATION"; ITER 
170  GOSUB  2990 
180  GOSUB  2330 
190  GOSUB  2190 

196  IF  A$<>"A"  THEN  200 

197  DELTA=FNRAD<DELTA) : G 
OSUB  10000 

200  GOSUB  1630 

205  CLS:  PR I NT "WORKING": 

PRINT" ITERATION ";ITER 

210  IF  F$="FAIL"  THEN  33 

00 

220  GOSUB  1370 

230  ITER=ITER+1 

240  IF  F$=“FAIL"  OR  ITER 

>IMAX  THEN  3300 

250  GOTO  200 

260  ’**  CONSTRAINT  TEST* 
* 

270  G$="PASS":  NC=0 

280  FOR  P=1  TO  NPLV 

290  IF  H<P)=0  THEN  445 

300  II=P  : GOSUB  1230 

310  FCR  N=1  TO  NL 

328  FCON=-l 

330  FOR  K=1  TO  3 

340  FOR  J=1  TO  3 

350  FCON=FCON+G  <  K , J  >  *E  <  N 

» J)*E<N,K) 

360  NEXT  J 

370  FCON=»FCON+S<K)*E<N,K 

> 

380  NEXT  K 

390  IF  FCON>0  THEN  G$="F 
AIL":  RETURN 

400  IF  FCON< -E5  THEN  440 


Comments 

20-40  commands  to  configure  the 
machine 


50  Implicit  integer  and  double 
precision 

80-90  convert  radians  to  degrees 
and  degrees  to  radians 

130  -  Gosub  input 

170  -  Gosub  invariants 

180  -  Gosub  transformations 

190  -  Gosub  initial  feasible  pt. 

200  -  Gosub  direction 

220  -  Gosub  new  thickness 


290  -  if  ply  thickness  zero, 
ignore  constraint 


300  -  Get  G  matrix  for  ply  being 
tested 


320  - 

G  .  .  e  .  e 
13  i 


380 

,  + 
3 


Solve 
G.  e  . 


l  l 


FCON  = 
1 


410  -  430  If  FCON  is  close  to 
zero  identify  constraint  as 
active,  make  list  in  C%  and  in¬ 
crement  constraint  counter 


97 


k  **  \ 


A  A.i  -V 


-V  -1  . 


A. 


410  NONC+1 
420  C5i<HC>  1  >=P 
430  C5i<HC>2)=N 
440  NEXT  N 
445  NEXT  P 
450  RETURN 
456  STOP 

460  ’**  GRADIENT  ** 

475  UNORM=0 

480  II=P:  GOSUB  1230 

490  FOR  L=1  TO  NPLS' 

500  IF  H<L)=0  THEN  700 

510  II=L:  GOSUB  1120 

520  FOR  J=1  TO  3 

530  R< J)=0 

540  FOR  K=1  TO  3 

550  R<J)=R(J)-Q<J,K)*E<N 

, K) 

560  NEXT  Kj  J 

570  FOR  J=1  TO  3 

580  V<J>=0 

590  FOR  K=1  TO  3 

600  V<J>=V(J>+AI<J,K)*R< 

K) 

610  NEXT  K» J 

620  2<L>=0 

630  FOR  J=1  TO  3 

640  FOR  K=l  TO  3 

650  2<L>=Z<L>+G<J,K>*<V< 

J)*E<NjK)+E(NjJ)*V<K>) 

660  NEXT  K 

670  2<L>-Z<L)+S< J)*V( J) 
680  NEXT  J 

690  UN0RM=UN0RM+2  <  L ) *2  <  L 

> 

700  NEXT  L 

710  UNORM=SQROJNORM) 

720  FOR  L=1  TO  NPLV 
730  Z<L)=2<L)/UNORM 
740  NEXT  L 
760  RETURN 
770  ’**  STRAINS  ** 

780  DIM  F<3,3) 

790  FOR  1=1  TO  3 

300  FOR  J=1  TO  3 

310  F<IjJ)=A<IjJ)+D<IjJ> 

*S 

820  NEXT  J,  I 

330  DET#=F< 1 j 1 )*F<2j  2>*F 
<3,3)+2*F<l,2)*F<2,3>*F< 
1,3)-F<2,2>*F(1,3)*F<1,3 
)-F<l j 1)*F<2j3)#F<2j3)-F 
<3j3)*F<1j2)*F(1j2) 

340  Aia.l)  =  <F<2.2)*F<3, 
3>-F<2j3)*F(2,3>)/'DET# 
350  AI<2,2)=<F<1,1>*F<3, 
3)-F<  1 » 3>*F< 1 j  3) >/DET# 
860  AI(1,2)=<F(1,3)*FC2j 
3)-F<l j2)*F<3j3))/DET# 
870  AI<3,3>=<F<1, 1>*F<2, 
2>-F<l,2)*F<l,2>>/DET# 
380  AI<1j3)=<F<1j2)*F<2j 
3>-F<2, 2>*F<  1 , 3>  VDET# 
390  AI<2»3>=(F<1 j2)*F<1 j 
3)-F<l.l)*F<2,3))/DET# 
900  AI<2,1>=AI<1,2>:AI<3 
,2)*AI<2,3):AI<3j1)=AI<1 


480 

ply 


-  Get  G  matrix  for  designated 


510  -  For  each  ply,  get  Q  matrix 
540 


560  R  =  -  TT-  A  e 

dn 


580  -  610  Y  =  | A-1 |  R 


Y  =  — —  E 
3hk 


620  - 

680 

$  (FCON)  =  fGi j ( ei 

+  -LeA 
9\ 

Ej) 

+  G.  (|£i) ]  Ak 

690  - 

730 

Normalize  V(FCON) 

790  -  820  "F"  is  the  A  matrix 
corresponding  to  a  point  S 
along  the  Z  vector 

830  - 

900 

invert  A 

920  - 

970 

Solve  t= 1 A  1 1  N  for 

each  independent  loading 


910  ERASE  F 

920  FOR  1*1  TO  NL 

930  FOR  >1  TO  3 

940  E<1» J>*0 

950  FOR  K*1  TO  3 

960  E<I. J)*E<r,J)+AI<J,K 

)*XN<  I  >  K> 

970  NEXT  K,J, I 

980  RETURN 

990  ’#*  A  MATRIX  ** 

1000  FOR  1=  1  TO  3 
1010  FOR  J=1  TO  3 
1020  A<I.J>=0:  D<I » J)=Q 
1030  NEXT  J,  I 
1040  FOR  1=1  TO  HPLV 
1050  1 1=1 •  GOSUB  1120 
1060  FOR  J=1  TO  3 
1070  FOR  K=1  TO  3 
1080  A<J>K)=A<J>K)+Q<J>K 
)#H< I ) 

1090  D< J>K)=tXJ>K)+Q< Jf K 

1100  NEXT  K,J,I 
1110  RETURN 
1120  '**  FORM  Q  ** 

1130  Q< 1 f 1 )=C< II > 1 ) 

1140  QC1»2)=C<II»3/'1 
1150  Q<1»3)=C<II»5) 

1170  Q<3> 1)=C<II»5) 

1180  Q<3»2)=C<II>6) 

1190  GK3,3>=C<II,4> 

1195  Q<2>3)=C<II>6) 

1200  Q<2,2)=C<II,2) 

1210  Q<2f 1)=C<II>3) 

1220  RETURN 
1230  ’**  FORM  G  ** 

1240  G<1> 1)=B<II> 1) 

1250  G<1»2)=B<II>3) 

1260  G< 1 > 3)=B< 1 1 1 5) 

1270  G<2»1)=B<II,3) 

1280  6<2>2)=B<II >2) 

1290  G<2»3)=B<II>6) 

1300  G<3» 1)=B<II>5) 

1310  G<3>2)=B<II»6) 

1320  G<3»3)=B<IIj4) 

1330  S<l)=B<IIi7) 

1340  S<2)=B<II»8) 

1350  S<3)=B<II,9) 

1360  RETURN 

1370  ’**  NEW  H  UECTOR  ** 
1380  SMAX=1E10 
1390  FOR  1=1  TO  NPLV 
1400  IF  Z<I)<>0  THEN  S=- 
H<I)/2<I> 

1410  IF  S>0  AND  S< SMAX  T 
HEN  SMAX=S 
1420  NEXT  I 
1430  F*="" 

1440  IF  SMAX>  10  THEN  F t 
=MFAIL"s  RETURN 
1450  S1=0:  S2=SMAX:  S=SM 
AX 

1460  IF  NC=0  THEN  1598 
1470  GOSUB  770:  GOSUB  26 
0 

1480  IF  G#=MFAIL"  THEN  S 

2=S  ELSE  S1=S 

1490  IF  SI -SMAX  THEN  153 

5 


1000  -  1100  The  matrix  D  is 
formed  so  that  along  the  Z  vector 

|A|  =  | A |  +  |D|  *  S 

where  S  is  a  scalor 


1130  -  1210  Convert  C  array  into 
3  x  3  Q  matrix  for  ply  designated 
by  II 


1240  -  1350  Convert  B  array  into 
3  x  3  G  matrix  for  ply  designated 
by  II.  Linear  failure  terms 
placed  in  vector  S 


1380  -  1420  Find  distance  along 
Z  to  find  h^  =  0  constraint 

1450  -  1500  Bisection  method  to 
find  distance  to  next  constraint. 
If  no  constraints  violated  at 
S  =  SMAX  then  stop  search 
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1500  S*<S1+S2V2 
1510  IF  S2-SKE2  AND  Sl  = 
0  THEN  F$="FAIL":  S=0:  G 
OTO  1650 

1520  IF  S1AS2-S1X4  THE 

N  1470 

1530  S=S/2 

1535  SREF=0 

1540  FOR  1=1  TO  NPLV 

1550  H<I>=H<I>+2<I>*S 

1560  IF  HCIXE2  THEN  H<I 

>=0 

1570  SREF=SREF+H<I)*H<I> 
1580  NEXT  I 

1590  S=0s  SREF=SQR ( SREF ) 

1600  0OSUB  990:  GOSUB  77 

0:  G0SU8  2020 

1610  IF  SREF-S<E2  THEN  F 

$="FAIL“ 

1620  FOR  1=1  TO  NPLV 
1630  H<I>=H<I>*S/'SREF 
1640  NEXT  I 
1650  S=0 

1660  GOSUB  990:  QOSUB  77 

0:  GOSUB  260 

1670  RETURN 

1680  ’**  DIRECTION  ** 

1690  2=0:  UNGRM= 1 

1700  FOR  1=1  TO  NPLV 

1710  X<I>=0 

1720  2=Z+SGN<H<I>) 

1730  NEXT  I 
1740  Z=XSQR<2> 

1750  IF  NC=0  THEN  1860 
1760  FOR  1=1  TO  NC 
1770  P=CX<I, l):  H=CX< I , 2 
) 

1780  GOSUB  460 

1790  FOR  J=1  TO  NPLV 

1800  LET  X<J>=X<J>-2<J> 

1810  NEXT  J,I 

1315  UNORM=0 

1820  FOR  J=1  TO  NPLV 

1830  UNORM=UNORM+X  <  J ) #X  < 

J) 

1840  NEXT  J 

1850  UNORM=SGR  <  UNORM ) :  T 
EST=0 

1860  FOR  1=1  TO  NPLV 
1870  X<I)=XaVUNORM 
1880  TEST=TEST+X< I >*2*SG 
N<H<I)) 

1890  NEXT  I 

1900  UNORM=0 

1910  FOR  1=1  TO  NPLV 

1920  2< I)=X< I >-TEST*2*SG 

N<H<I>) 

1930  UNORM=UNORM+2 < I )*2< 
I) 

1940  NEXT  I 

1950  IF  UNORM< IE-6  THEN 
F*=“FAIL":  RETURN  ELSE 
F«=“  " 

1960  UNORM=SQR< UNORM) 
1970  FOR  1=1  TO  NPLV 
1980  2<I)=2<I>/UN0RM 
1990  NEXT  I 
2000  GOSUB  990 
2010  RETURN 


1530  -  1600  at  point  halfway 
between  constraints,  use  strain 
ratio  routine  to  find  how  much 
the  laminate  thickness  can  be 
reduced 


1610  If  change  in  thickness 
small,  set  flag  to  halt  program 


1620  -  1660  Update  h  vector,  A 
matrix,  strains 


1760  -  1840  For  each  active 
constraint  call  gradient  sub¬ 
routine.  Sum  negative  of  each^ 
gradient  into  X  and  normalize  X 


1860  -  1890  Take  dot  product  of 
X  and  unit  normal  to  Ehi  =  const, 
plane 


1910  -  1940  2  is  a  vector  parallel 

to  the  Ehi  =  const,  plane  and 
pointing  away  from  the  active 
constraints 


1950  if  the  magnitude  of  Z 
is  very  small,  a  local  minima  has 
been  reached 


V, 

y 


$ 


2020  '**  STRAIN  RATIO  ** 
2030  FOR  P=1  TO  NPLV 
2040  IF  H<P>=0  THEN  2160 
2050  II=P:  GOSUB  1230 
2060  FOR  N=1  TO  NL 
2070  B#=0:C#=0 
2080  FOR  1=1  TO  3 
2090  FOR  J=1  TO  3 
2100  C#=C#-SREF*SREF*Gd 
, J)*E<N,I)*E<N,J) 

2110  NEXT  J 

2120  B#=B#-SREF*Sd)*E<N 

,1) 

2130  NEXT  I 

2140  SUAL=<-B#+SQR  <B#*B 
#-4*C#*< 1-E6>  >  >A2*d-E6 
)) 

2150  IF  SUAL >S  THEN  S=SU 
AL 

2155  NEXT  N 
2160  NEXT  P 
2180  RETURN 
2190  ’**  IFP  ** 

2200  Z=1''SQR<NPLV) 

2210  FOR  1=1  TO  NPLV 
2220  Zd)=Z:  H< I )=Z 
2230  NEXT  I 
2240  GOSUB  999 
2250  S=0:  SREF=1 
2260  GOSUB  770:  GOSUB  20 
20 

2270  FOR  1=1  TO  NPLV 
2280  H< I )=H< I )#S 
2290  NEXT  I 
2300  S=0 

2310  GOSUB  990:  GOSUB  77 

0:  GOSUB  260 

2320  RETURN 

2330  '**  TRANSFORM  ** 

2340  FOR  1=1  TO  NPLV 

2350  C2=C0S<2*Td>>:  C4= 

C0S<4*Td>> 

2360  S2=SIN<2*Td)):  S4= 
SIN<4*T<I)) 

2370  Bd,  l)=Ud>+C2*U<2) 
+C4*U<3> 

2380  B<I,2)=U<1)-C2*U<2) 
+C4*U<3> 

2390  Bd >  3)=U(4)-C4#U<3> 
2400  Bd»4)=U<5)-C4*U<3) 
2410  Bd,5)=S2/'2*U<2)+S4 
#U<3> 

2420  Bd,6)=S2/'2*U<2>-S4 
*U<3> 

2430  Bd,7)=U<6>+C2*U<7> 
2440  Bd,8>=U<6>-C2*U<7> 
2450  Bd,9)=S2*U<7> 

2460  Cd,  1 >=Ud  >+C2*U<2> 
♦C4*U<3> 

2470  Cd  ,2>=U< 1 )-C2*U<2> 
+C4*U<3> 

2480  Cd,3>=U<4>-C4*U<3> 
2490  Cd >  4)=U<5)-C4*U<3> 
2500  Cd > 5)=S2-/2*U<2)+S4 
*U<3> 

2510  C<I»6)=S2/2*U<2)-S4 
*U<3> 

2520  NEXT  I 
2530  RETURN 


2030  -  2140  For  each  possible 
constraint  solve  for  S  in 

G..C.S.  iSREFif.  +  G  £  .(SREF), 
13  l  3  „2  ii  S 


-  1  =  -E6 


2150  Take  smallest  value 
(corresponds  to  closest  constraint) 


2200  -  2310  For  equal  ply  ratios, 
find  the  smallest  laminate  thick¬ 
ness  which  does  not  violate  any 
constraints.  Initialize  A  matrix, 
strains,  and  constraint  list 


2370  -  2450 

Transform  failure 

parameters  : 

in  following  order 

b(i,i)«gu 

B(I,5)=G16 

B(I,2)=G22 

B(I,6)=G26 

B(I,3)=G12 

B(I,7)=G1 

b<i'4>°666 

B(I,8)=G2 

B(I,9)=G3 

2460  -  2510 

Transform  modulus 

in  following  order 

C  ( 1 , 1 )  =  Qn  C  ( 1 , 5)  =  Q16 
C(I,L)=Q22  C(I,6)=Q26 

C(I,3)=Q12 
C  ( 1 , 4 ) =Q, , 
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2540  ’***  INPUT  *** 

2550  CLS 

2600  PRINT  "PRESS  ANY  KE 
V  WHEN" SPRINT  "DESIRED  M 
ATER I AL “ : PR I NT " APP  EARS " 
2610  FOR  K=1  TO  750: NEXT 
2620  FOR  M=0  TO  6 
2640  IF  M=6  THEN  Mf="NEU 
MATERIAL"  ELSE  GETXM.EX 
, EY.UX.ES.TPLYjXT  > VT , xc . 
YC.SS.M* 

2650  CLS: PRINT  M*: SOUND 

20,1 

2660  FOR  J=1  TO  200 
2670  IF  INKEYSO""  THEN 
2700 

2675  NEXT  J,M 

2630  GOTO  2620 

2700  IF  M=6  THEN  GOSUB  9 

000: GOTO  2600 

2705  CLS  .‘PRINT  "f  ":M#  J" 

1" 

2710  PRINT  "HOW  MANY" 
2720  INPUT  "PLY  GROUPS": 
NPLV 

2730  CLS:  PRINT  "ENTER  P 
LY  GROUP" 

2740  PRINT  "ORIENTATIONS 

II 

2750  FOR  1=1  TO  200 
2760  NEXT  I 
2770  CLS 

2780  FOR  1=1  TO  NPLV 
2790  PRINT  "PLY  ”5  I 
2800  INPUT  T< I > 

2810  T < I )=FNRAD<T < I) ) 
2820  NEXT  I 

2830  PRINT  "ENTER  NUMBER 
OF" 

2840  PRINT  "INDEPENDENT 
LOAD" 

2850  INPUT  “CONDITIONS"; 
NL 

2900  FOR  1=1  TO  NL 
2910  CLS: FR I NT  "LOAD  ":I 
;  "  IN  MPa.  “ 

2920  INPUT  "Nl=";XNd,l> 
2930  INPUT  "N2=";XN<I»2> 
2940  INPUT  "N6="*XN(I»3> 
2950  FOR  J=1  TO  3 
2960  XN<I, J)=XN< I>  J)*1E6 
2970  NEXT  J, I 
2980  RETURN 

2990  ’**  INUARIENTS  ** 
3050  UV=1/<1-UX*UX*EY/EX 
) 

3060  GXX=UY*EX*1E9:  QYY= 
UV*EV*1E9 

3070  QXV=UV*UX*EY*1E9:  Q 
S=ES*1E9 

3080  U<i>=<3*QXX+3*QVV+2 
*QXV+4*QS)/'B 
3090  U<2>=<QXX-QYYV2 
3100  U  <  3  > »  <  GXX+Q YV  -2*QXY 

-4*QSV8 


2600  -  2675  List  available 
materials. Get%  is  an  HX-20 
command  to  get  data  from  a  non 
violatile  RAM  file 


3050  -  3280  Calculate  invariants 
for  use  in  transformations.  Note 
that  some  variables  like  EX  and 
EY  get  reused,  so  their  value  may 
not  be  what  you  might  expect  after 
routine  is  called 


3110  U<4)*<QXX+QVV+6*QXY 
-4*QS)/8 

3120  U  <  5  > = <  QXX+GYV-2*QXV 
+4*QS)/S 

3130  EX=1E-12'/<XT*XC>:  E 
V=1E-12/'<VT*VCV.  ES=1E-1 
2/<SS*SS) 

3140  FX=<1/XT-1/XCV1E6: 

FV=<1/VT-1/'VCV1E6 
3150  EXV=-SQR  <EX*EY>/2 
3160  GXX=EX*QXX*QXX+2*EX 
Y*QXX*0XV+EY*QXY*QXY 
3170  GVV=EX*G!XV*QXY+2*EX 
v*qxv*ciyv+ev*qvy*qyv 
3180  GXY=EX*0XX*G!XV+EXV* 

<  QXX*Q Y V+QXV  *QXV >+EY*QXY 
♦QVY 

3190  GSS=ES*QS*QS 
3200  GX=FX*QXX+FV*QXV 
3210  GY=FX*QXV+FV*QYV 
3220  ua>=<3*GXX+3*GYY+2 
*GXY+4*GSS>/'8 
3230  U<2)=<GXX-6VYV2 
3240  U<3>=<6XX+GVV-2*GXY 
-4*GSSK3 

3250  U<4)-<GXX+GYV+6*6XV 
-4*GSS)/8 

3260  U  <  5  ) = <  GXX+GVV-2*GXV 
+4*GSS)/8 

3270  U(6)=<GX+6V)/2 
3280  U<?)=<GX-GYV2 
3290  RETURN 
3300  ’**  OUTPUT  ** 

3302  SOUND  15, 2:SOUND50, 
2 

3305  K*="Hit  any  key  to 
cont":U*="MN/fn" 

3310  CLS:  TEST=0 
3320  FOR  1=1  TO  HPLV 
3330  TEST=TEST+H<I>s  NEX 
T  I 

3350  PRINT  "TOTAL  THICKN 
ESS=" 

3360  PRINT  TEST;  m  m.  " 
3370  PRINT  USING  "####.  i 
*  Plies"; TEST /TPLV 
3375  PRINT  K»; 

3380  A*=INKEY*:IF  A$<>"" 
THEN  3380 

3390  IF  INKEV**""  THEN  3 
390 

3400  CLS:PRINTMPress  V  i 
f  printout", "of  displays 
d  result  is  desired.  Pr 
ess  Nif  not"; 

3410  FOR  1=1  TO  600: NEXT 
I 

3415  Af=INKEY*:IF  A*<>" " 
THEN  3415 

3420  CLS: RESTORE  6120 
3425  JJ=0:A*= INKEY* 

3430  FOR  1=1  TO  8 
3440  READ  A*: CLS: PRINT: P 
PINT  A*: SOUND  20,1 
3443  A*=INKEV*: IF  A*="“ 
THEN  3445 

3450  PRINT  A*;:FOR  KK=1 
TO  75: NEXT  KK 


3455  IF  A$“"VH  THEN  JJ=J 
J+i:C*'.<JJ,l>  =  I 
3460  NEXT  I 

3464  IF  JJO0  THEN  LPRIN 
T  STRING*<24, "i" > 

3465  FOR  KK=1  TO  JJ 
3470  ON  CX<KK» 1 )  GOSUB  5 
000,5200. 4000,4200,4400, 
4600,4800,7500 

3485  LPRINT 
3490  NEXT  KK 

3495  CLS : PR I NT " F IN I SHED " 

3496  IF  INK£VS=" “THEN349 
6  ELSE  RUN 

4000  ’**  PLV  RATIO** 

4002  CLS: LPRINT  "Total  t 
hickness=“ 

4004  LPRINT  USING  ".#### 

AAAA  m  « • Jggl 

4006  LPRINT  USING  ”####. 
##  Plies"; TEST /TPLV 
4008  LPRINT 

4030  A*=" ANGLE  RATIO  # 
PLIES" 

4040  LOCATE  0,1:  PR I NT  A* 
: LPRINT  A* 

4050  FOR  1=1  TO  NPLV 
4060  A=CINT<  <FNDEG<T<  D ) 
*1E2))/1E2 

4070  B=CINT<<H<I)/TEST*1 
E4))/1E4 

4080  C=C I NT ( <  H ( I V  TPL V* 1 
E2))/1E2 

4090  PRINT  A;TAB<6);B;TA 
B<13)!C 

4100  LPRINT  a;tab<6);b;T 

AB<13>5C 

4120  NEXT  I 

4150  RETURN 

4200  ’**  STRENGTH** 

4210  LPRINT  "STRENGTH  RA 
TIOS" 

4215  LPRINT  H1=ULTIMATE 
CTPOTKI”  ! 

4220  LPRINT  ">1  IS  SAFE" 
4225  FOR  1=1  TO  NL 
4230  LPRINT  "LOADING  "I 
4235  LPRINT  "PLV" 

4240  FOR  P=1  TO  NPLV 
4245  IF  H<P>=0  THEN  4305 
4250  1 1 =P: GOSUB  1230 
4255  A#=0: B#=9 
4260  FOR  J=1  TO  3 
4265  FOR  K=i  TO  3 
4270  A#=A#+G<J, K)*E< I ,  J) 
*E<I,K> 

4275  NEXT  K 

4280  B#=B#+S< J)*E< I, J) 

4235  NEXT  J 

4290  A#=<-B#+SQR<B#*B#+4 
*A*))/<2*A#> 

4295  A=FIX<A#*1E4+.  5)/lE 
4 

4300  LPRINT  FNDEGCT'P) ) ; 
TAB(10)*A 
4305  NEXT  P, I 
4310  RETURN 


3400-3490  Branch  tor  various 
output  routines 


4000-7560  Output  routines  and 
laminate  anal ys i s 


4210  -  4305  Strength  ratio  is 
defined  as  the  value  of  R  in 

G.  .^E.R2  +  G.C.R  -1  =0 
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4400  ’**STRAINS** 

4410  LPRINT  TAB<4) 5 "LAM I 
HATE  STRAINS” 

4420  FOR  N=1  TO  NL 
4430  LPRINT  "LOADING  "N 
4440  LPRINT  USING  "el=+# 
.###E-03";E<N,1>*1£3 
4450  LPRINT  USING  "e2=+# 
. ###E-03" >E<N>2)*1E3 
4460  LPRINT  USING  "e6=+# 
.  ###E- 03 " *  E  <  N  >  3  >  * 1 E3 
4465  NEXT  N 
4470  RETURN 
4600  ’  **A  MATRIX** 

4610  CLS 

4620  LPRINT "Norm.  |A|  in 
GPa.  " 

4630  FOR  1=1  TO  3 
4640  FOR  J=1  TO  3 
4650  D<I>J)=A<I>  J)-'1E9'T 
EST 

4660  NEXT  J, I 
4670  GOSUB  7000 
4680  RETURN 
4800  ’A  INUERSE 
4810  L PR I NT “Compliance 
(normal i red) " 

4820  LPRINT" in  1/TPa. " 
4830  FOR  1=1  TO  3 
4840  FOR  J=1  TO  3 
4850  D<I>J)=AI(I> J)*TEST 
*1E12 

4860  NEXT  J,I 
4870  GOSUB  7000 
4830  RETURN 

5000  LPRINT  "Ma+er ial  Pr 

npor  +  i  pc 11 

5010  GET*M,EX,EV>UX,ES,T 
PLVj  XT  j  VT » XC j  VC>  SS>  M$ 
5015  LPRINT  Mi 
5020  LPRINT  "EX=";EX; "GP 
a" 

5030  LPRINT  "EV=";EV;"GP 
a" 

5040  LPRINT  "ES=";ES;"GP 
3" 

5050  LPRINT  "UX=" >  UX 
5060  LPRINT  "X="!XT* "MPa 

H 

5070  LPRINT  "X’="5XC;"MP 
a" 

5072  LPRINT  "V=";VT;“MPa 

II 

5074  LPRINT  "V’ =" 5  VC; "MP 
a" 

5030  LPRINT  "S=";SS;"MPa 

II 

5090  LPRINT  "Ply  Thickne 
ss">  TPLV'm" 

5095  RETURN 
5200  ’ **LOADS** 

5210  FOR  1=1  TO  NL 
5220  LPRINT  "LOADING  "?I 
5230  FOR  J=1  TO  3 
5240  Ai=STR*(J): IFJ=3THE 
N  A$="  6" 

5250  LPRINT  ”N"*A$; "="5X 
NCI, J)/lE6;Ui 
5260  NEXT  J.  I 
5270  RETURN 


6000  DATA  10.5E-5?.  1>1E- 
6 

6120  DATA  Ply  properties 
, Loads » Total  thickness  & 
Ply  rati  os, Strength 
ratios 

6130  DATA  Laminate  strai 
ns? Stiffness  matrix? Comp 
liance  matrix? Engineer  in 
9  const. 

6500  ’**  ANAL VS IS  ** 

6505  CLS 

6510  TEST*0 

6520  FOR  1=1  TO  NPLV 

6530  PRINT  "NUMBER  OF  PL 

IES  AT"?FNDE6<T<I>  > ? "DEQ 

PEES"!  INPUT  H<I> 

6540  H< I >*H< I >*TPLV: TEST 
=TE3T+H<I) 

6550  NEXT 
6555  CLS 

6560  3=0:GOSUB  2990: GOSU 
B  2330: GOSUB  990:GOSUB  7 
70 

6570  GOTO  3300 
7000  ’ FANCV 

7010  LPRINT  ", - 1 — 


7020  LPRINT  USING  "I###. 
###">D< 1? 1 >?  D<1?2>  ?D<1?3 


7040  LPRINT  AS 
7050  LPRINT  USING  "I###. 
###";D<2.  l).D<2,2)?0<2?3 
) 

7060  LPRINT  AS 
7070  LPRINT  USING  "I###. 
###"?D<3? 1>?D<3?2)»D<3?3 
> 

7080  LPRINT  - 1 — 

_ i - 1  " 

7100  RETURN 

7500  ’**  ENG.  CONST.  ** 
7505  LPRINT  "ENGINEERING 
CONSTANTS": LPRINT 
7510  LPRINT  USING" El=### 
.  #  GPa"?  l^AI  <  1 » 1  )/TEST^'l 
E9 

7515  LPRINT  USING"E2=### 
.#  GPa";lxAI<2.2)/TESTxl 
E9 

7520  LPRINT  USING"E6=### 
.  #  GPa":l/AI<3?3)/TEST/l 
E9 

7525  LPRINT  USING"E2=### 
.  #  GPa"!  lxAI  <2?  2)/TEST?"l 
E9 

7530  LPRINT  USING"v21=## 
.###"; -AI<1,2VAI<1?1> 
7540  LPRINT  USING"v61Mt 
.  ###":Aia?3VAI<l?l) 
7550  LPRINT  USING" o 16=## 
. ###"*  AI< 1 »3)^AI(3>3) 
7560  RETURN 


•6500-6570  MAIN  tor  perform  in 
only  laminate  analysis  (numbe 
of  plies  is  a  given) 


I  Ob 


9000  PRINT"REUIEW  OR  NEW 

II 


9010  INPUT  "DATA  <R/N>"; 

A* 

9020  IF  A*="R"  THEN  9190 
9030  PRINT "WHICH  MATERIA 
L":PRINT"WILL  VOU" 

9040  INPUT  "REPLACE  <0-5 
)"JI 

9050  INPUT  ” EX < GF'a >  =  "  : EX 
9060  INPUT  "EV<QPa>="*EV 

9079  INPUT  "UX=";UX 
9075  INPUT  "ES<GPa>=";ES 

9080  INPUT  "X< MPa )=";X 
9090  INPUT  "X?  <MPa)=" ; XX 
9100  INPUT  hY<HP3>=";V 
9110  INPUT  "V’ <MPa>="5  VV 
9120  INPUT  "S<MPa)=" jS 
9130  INPUT  "THICKNESS  <m 
)=" ’ TPLY 

9140  INPUT  "NAME  <15  CHR 
.  MAX.  ) " ; M* 

9150  PUT*4I , EX.  EV, UX, ES,  T 
PLV.X, Y.XX.VV.S.Mf 
9160  PRINT  "ADDITIONAL": 
INPUT  "CHANGES  < WN> " 5 A$ 
9170  IF  A$="V"  THEN  9000 
9180  RETURN 

9190  PRINT"REUIEW  WHICH" 
: INPUT "MATERIAL  <8-5)"?M 
9200  GOSUB5000 
9210  GOTO  9160 
9500  OPEN  “ I ” j  # 1 . " CAS 1 : D 
AT  A" 

9510  FOR  1=0  TO  5 

9520  INPUT  #1 > EX, EV>  UX, E 

S,T,X.Y,XX,YY,S,M* 

9530  PUTX I , EX. EV, UX . ES ,  T 


,  X.  V,  XX,  VY.  3,  Mt- 
9540  NEXT 
9550  CLOSE  #1 
9560  DELETE  45 
9570  GOTO  50 
10000  'THETA 
10010  L=0 
10020  SPEF=0 
10030  FOR  1=1  TO  NPLV 
10040  SREF=SREF+H< I >*H< I 
? 

10050  NEXT 

10060  FOR  P=1  TO  NPLV 
10070  C2=2*C0S<2*T<P>> 
10080  C4=4*C0S<4*T <P> > 
10090  S2=2*S I N < 2*T  < P  >  > 
10100  S4=4*SIN<4*T<P>> 
10110  D<1. 1 >=-U<2)*S2-U< 
3>*34 

10120  D<2.2)='J<2)*S2-U<3 
)*S4 

10130  D<3, 3>=U<3)*S4 
10140  D<  1 . 3)  =U<2)*C2,/2+U 
<3>*C4 

10150  D<2.3)=U<2)*C2/2-U 
<3>*C4 


9000-921 0 
properties 


Enter  material 
into  library 


9500—9570  Routine  -for 
automatically  entering  material 
properties  from  a  cassette 
tape.  To  use,  line  45  should 
read  GOTO  9500,  and  load 
prociram.  With  tape  still 
connected,  run  program  and  the 
properties  will  load. 


10000  Angle  optimization 
subrou  tine 


101 10-10190 
derivatives 
pajama  ters 


An  gu 1 ar 
o+  -failure 


107 


10160  D(1,2)=D<3,3>:D<2, 
l)=D(l>2)!D<3>2)=D(2>3)s 


D<3> 1  )=D< 1>3) 

10170  X<1>=-U(7>*S2 
10180  X<2)=-X<1) 

10190  X(3>=U(?>*C2 
10200  I I=P: 60SUB  1230 
10210  FOR  N=1  TO  NL 
10220  FCON=-l 
10230  FOR  K=1  TO  3 
10240  FOR  J=1  TO  3 
10250  FCON=FCON+G<J,IO*E 
(N,J)*E<N,IO 
10260  NEXT  J 

1 0270  FCON*FCON+S <  K > *E (  N 
>  K) 

1O280  NEXT  K 
10290  L=L+1 
10300  CON(L)=FCON 
10310  FOR  PF-1  TO  NPLV 
1 0320  C2=2*C0S ( 2*T < PP  >  > : 
C4=4*C0S(4*T<PP>> 

1 0330  S2=2*S I N  ( 2*T < PP >  > : 
S4=4*SIN<4*T <PP>) 

1 0340  h  <  1 , 1  )  =-IJ  <  2  >  +82-U < 
3)  *84 

10350  A<2> 2>=U(2>*S2-lK3 
)*S4 

10360  A<3>  3>=U<3)#S4 
1 0370  A  < 1 , 3 ) =U  < 2 > *C2-2+U 
(3)*C4 

10380  A(2,3)=U<2>*C2'2-U 
C3)*C4 


10390  A< 1 » 2)=A<3»3) : A<2> 
1>=A<1,2>:A<3,2>=A<2,3>: 
A<3> 1  )=A( 1 >3) 


104O0  FOR  J=1  TO  3 
10410  R< J)=0 


10420  FOR  K=1  TO  3 
1 0430  R <  J  >  =R (J ) + A  <  J , K >  *E 
<N,K)*H<P) 


10440  NEXT  K,J 
10450  FOR  J=1  TO  3 
10460  V(J)=0 
10470  FOR  K=1  TO  3 


10480  V(J>=V<J)-AI < J»K>* 
R(K) 


10490  NEXT  K>  J 
10500  DUM=0 
10510  FOR  J=1  TO  3 
10520  FOR  K=1  TO  3 
10530  DUM=DUM+G<  J,  K >*<V< 
J  >  *E ( N  >  K) +E  <  N  >  J ) #V <  K  >  > 
10540  IF  P=PP  THEN  DUM=D 
IJM+D  (J»K)*E(N>  J)  *E  <  N ,  K  > 
10550  NEXT  K 


10560  DUM=DU«+S<J)*V<J> 
10570  IF  P=PP  THEN  DUM=D 
UM+X  ( J ) *E  <  N  >  J ) 

10580  NEXT  J 
10590  UK  L  >  PP  >  =DUM 
10600  NEXT  PP 
10610  NEXT  N>  P 
10620  ZHAX=0 
10630  FOR  P=1  TO  NPLV 
1 0640  DUM=0 : 0UM2=8 : DUM3= 
0 


10250  10270  Calculate  va 
failure  equation  for  each 


10340-10390  An  gu 1 ar 
derivatives  of  Q  matrix 
t  e  rms 


1 040 0  —  1 0490  So  1 v e  f  or 
partials  of  strain  with 
respect  to  angle 


10510—10600  Solve 
equation  (29) 


10  630  10  700  S  o  v  e  for 
equation  (28) 


u  e  o  f 
1  oad 


1U8 


10650  FOR  J=1  TO  L 
1 0660  DUM=DUM+CON < J >  *W ( J 
,P) 

10670  DUM2=BUN2+W<J,P> 

1 O680  WJM3=DUM3+C0N < J > 
10690  NEXT  J 

1 0700  ZCP> =DUM-DUM2*DUM3 
/NPLV 

10710  IF  ABS<Z<P))>ZMAX 
AND  Z<P)<>0  THEN  ZMAX=AB 
S<Z<P>> 

10720  X<P>=0 
10730  NEXT  P 
1O740  FOR  1=1  TO  NPLV 
10750  IF  ZMAX=0  THEN  RET 
URN  ELSE  Z<I>=-Z<I VZMAX 
10760  NEXT  I 
10770  T=l:TEST=SREF 
10780  FOR  1=1  TO  NPLV 
1 0790  CON ( I > =X < I ) : X ( I > =Z 
a  >*T 

10300  YX I )=CINT (X< I ) > 
10310  T< I ) =T ( I >  +  < X < I > -CO 
N( I ) >*DELTA 
18320  NEXT 

10830  GOSUB  2330: 3=0: GOS 
UB  990: GOSUB  770 
10840  GOSUB  2020 
18850  IF  S<TEST  THEN  T=T 
+  1 : TEST=S: GOTO  18780 
10860  FOR  1=1  TO  NPLV 
10370  T  < I >  =T (I > - < X < I ) -CO 
N( I > )+DELTA 
10330  NEXT 

18390  S=0: GOSUB  2330:603 
UB  990: GOSUB  770: GOSUB  2 
820 

10900  FOR  1=1  TO  NPLV 
18910  H  <  I  >  =H  ( I  >  *S/SREF 
10930  1 1 =P: GOSUB  1230 
10940  NEXT 

10960  IF  T=1  THEN  RETURN 
ELSE  GOTO  10000 


10740-10760  Normalize  Z 
by  largest  component 


107S0 -10820  Incremental 
step  of  all  angles 


10830  Update 
t  r an sforma t i on s  ,  strains, 
and  scale  total  thickness 


10860-10880  After  minimum 
past,  go  back  one  step 


1 0900-1 0940  Update  ply 
group  thickness- 


10  960  If  any  progress- 
made,  go  back  and  try  a 
new  direction.  If  not, 
return  to  ma i n 


1U9 


APPENDIX  C 


Orthotropic  Laminate  Optimization  Program 

The  following  program  produces  a  thickness  optimized  laminate  that 
is  constraint  to  be  orthotropic.  A  search  can  be  made  for  the  best 
orientation  of  the  orthotropic  axis.  In  the  final  result,  ply  angles 
are  measured  from  one  of  the  orthotropic  axes  (the  original  1  axis)  plus 
a  rigid  body  rotation  is  given.  Angles  appear  to  stay  constant,  but  the 
rigid  body  rotation  must  be  added  to  get  the  angle  to  the  laminate  1 
axis  (see  Figure  15).  The  failure  theory  is  based  on  the  strain  sphere 
approximation  of  the  first-ply  failure  inner-envelope.  The  laminate 
must  remain  balanced.  Instead  of  entering  both  the  plus  and  minus 
angle,  only  one  is  entered  and  the  program  assumes  the  presence  of  the 
other.  The  final  thickness  must  be  divided  evenly  between  a  plus  theta 
and  a  minus  theta  ply  group. 

Running  the  program  is  similar  to  running  the  program  listed  in 
Appendix  B.  The  only  differences  are  that  if  optimimum  orientation  is 
desired,  the  search  limits  and  maximum  error  must  be  entered.  The 
search  limits  are  the  angles  between  which  the  best  laminate  is  thought 
to  lie.  If  a  minimum  is  not  found  between  the  given  limits,  the  program 
automatically  extends  the  limits,  but  this  is  time  consuming. 

All  angles  (and  the  error)  are  entered  as  degrees. 


10  ’**  MAIN  CLASS** 

20  WIDTH  20, 10 

30  CLEAR  75,330 

40  DEFFIL  55,0 

50  DIM  XN<3> 3>»Q11<4),Q2 

2<4) , Q12<4) , Q66<4) , H<4) » 

R<3>,T<4> 

60  DIM  X<4),V<3),2c:4),E< 
3,3),Cf<4),U<5),U<5> 

70  DEF  FNRAD<X>=X/ 180*3. 
14159 

80  DEF  FNDEG<X)=X* 180/3. 
14159 

90  RESTORE 
100  READ  E2,E5,E6 
110  ITER=1 
120  G0SU8  1830 
138  GOSUB  2180 
140  GOSUB  1720 
150  INPUT  "OPT.  ROTATION 
<V/N)"5A$ 

160  IF  A$="V"  THEN  GOTO 
4180 

170  GOSUB  1580 

180  GOSUB  1160 

190  CLS: PRINT  "WORKING, I 

TERATION-JITER 

200  IF  F$="FAIL"  THEN  25 

0 

210  GOSUB  840 

220  ITER=ITER+1 

230  IF  F$= “FAIL"  THEN  25 

0 

240  GOTO  180 
250  H=0 

260  FOR  1=1  TO  NPLV 
270  H=H+H< I ) : NEXT  I 
280  GOTO  2330 
290  ’**  CONSTRAINT  TEST* 
* 

300  G*="PASS":  NC=0 

310  FOR  N=1  TO  NL 

320  FCON=CE<N, 1)*E<N, 1)+ 

E<N,2)*E<N.2)+E<N,3)*E<N 

,3)/2)/EMAX-l 

330  IF  FCOH>0  THEN  G$="F 

AIL": RETURN 

340  IF  FC0N<-E5  THEN  GOT 
0  370 

350  NC=NC+1 

360  C*<NC)=CHR*<N> 

370  NEXT  N 
380  RETURN 
400  ’**  GRADIENT  ** 

410  UNGRM=0 
420  FOR  L=1  TO  NPLV 
430  IF  H<L)=0  THEN  Z<L>= 
0JGOTO  520 

440  R< 1 >=-Ql 1<L)*E<N,1>- 
Q12<L>*E<N,2> 

450  R(2>=-G12<L)*E<N, 1 )- 
Q22<L)*E<N,  2) 

460  R<3>=-Q66<L)*E<N,3> 
470  V< 1 )=AI 1 1*R< 1 >+AI 12* 
R<2) 

480  V<2)“AI 12*R< 1 )+AI22* 
R<2> 


100  Error  and  numerical  offset 
cons tan ts 

120  Input 

130  Invar i ents 

140  Transformations 

160  Branch  for  optimum 
orientation 

170  Initial  feasible  point 

180  Direction 

210  New  position  in  design 
space 


300-380  Equation  (?) 


400-520  Partial  derivatives  of 
strain 


490  V<3>*AI66*R<3> 

500  Z<L)=2*V(l)*E(N,l)+2 
*V<2)*E<N»2)+V<3>#£(N>3> 

519  UNGRM=UNORM+Z<  L ) *2  <L 
) 

520  NEXT  L 

530  UNORM=SQR  <  UNORM  > 

540  FOR  L=1  TO  NPLY 
550  2<L>=2<L)/'UN0RM 
560  NEXT  L 
570  RETURN 
580  ’**  STRAINS  ** 

590  F11=A11+D11*S 

600  F12=A12+D12*S 

610  F22=A22+D22*S 

620  F66=A66+D66*S 

630  DET=F11*F22-F12*F12 

640  Alt  1=F22/DET 

650  AI22=F1 1/DET 

660  AI12=-F12/DET 

670  AI66=1/F66 

680  FOR  1=1  TO  NL 

690  E<I,l)=AIll*XNa,l>  + 

AI12#XN<I>2) 

700  E< I >2)=AI 12*XN< I j 1)+ 
AI22*XN< I >  2> 

710  E<I»3)=AI66*XNU,3> 

720  NEXT  I 

730  RETURN 

740  ’**  A  MATRIX  *+ 

750  A11=0:A22=0JA12=0.‘A6 
6=0 

760  Dll=0:D22=0:D12=0:D6 
6=0 

770  FOR  1=1  TO  HPLV 
730  AU=A11+Q11<I)*H<I>: 
Dll=Dll+Qll<n*Z<I) 

790  A22=A22+Q22C I ) *H< I ) : 
D22=D22+Q22< I >*Z< I } 

800  A12=A12+Q12< I >*H< I >: 

D12=D12+Q12(I)*Z<n 

810  A66=A66  +Q66< I )*H< I ) 

! D66=D66+Q66< I )*Z< I) 

320  NEXT  I 
830  RETURN 

840  REM  ***NEU  POSITION 
*** 

850  SMAX=1E10 

860  FOR  1=1  TO  NPLV 

870  IF  Z<I)<>0  THEN  S=-H 

<I)/Z<I) 

880  IF  S>0  AND  S<SMAX  TH 
EN  SMAX=S 
390  NEXT  I 
900  F$=“" 

910  IF  SMAX>  10  THEN  F$= 
"FAIL":  RETURN 
920  Sl=0s  S2=SMAX:  S=SMA 
X 

930  IF  HC=0  THEN  1070 
940  GOSUB  580:  GOSUB  290 
950  IF  G$="FAIL"  THEN  S2 
=S  ELSE  S1=S 

960  IF  S1=SMAX  THEN  1010 
970  S=<Sl+S2)-'2 
980  IF  S2-SKE2  AND  S1=0 
THEN  F*="FAIL":  S=0:  GO 
TO  1130 

990  IF  S1''<S2-S1  )<4  THEN 
940 

1000  S«S/2 


530-560  Normalize  oradient 


590-620  Update  A  matrix  -for 
point  S 

600-670  Invert  A  assuming 
ortho tropic  laminate 

680-710  Solve  tor  strains 


750-800  Form  A  matrix 

850-1000  Bisection  search  -for 
next  constraint 


1010-1140  At  haltway  point, 
rescale  1  ami  ate  and  update 
thickness  uec  tor- 


1010  SREF-0 
1020  FOR  1=1  TO  NPLV 
1030  H<I)=H<I)+2<IV*S 
1040  IF  HCIXE2  THEN  HCI 
)=0 

1 058  SREF=SR£F+H< I)WI) 
1060  NEXT  I 

1070  S=0:SREF=SQR(SREF> 

1080  GOSUB  740:  GOSUB  53 

0:  GOSUB  1510 

1090  IF  SREF-SCE2  THEN  F 

*="FAIL" 

1100  FOR  1=1  TO  NPLV 
1110  H<I>=H<I>*S/SREF 
1120  NEXT  I 
1130  S=0 

1140  GOSUB  740:  GOSUB  58 

0:  GOSUB  290 

1150  RETURN 

1160  ’**  DIRECTION  ** 

1170  W=0:  UN0RM= 1 

1180  FOR  1=1  TO  NPLV 

1190  X(I)=0 

1200  U=W+S6N<H< I ) ) 

1210  NEXT  I 
1220  W=1/SQR<W) 

1230  IF  NC=0  THEN  1350 
1240  FOR  1=1  TO  NC 
1250  N=ASC<C$ <I> ) 

1260  GOSUB  400 

1270  FOR  J=1  TO  NPLV 

1280  LET  X<J)=X<J)-Z<J) 

1290  NEXT  J,  I 

1300  UNORM-0 

1310  FOR  J=1  TO  NPLV 

1320  UNORM=UNORM+X  <  J > *X  < 

J) 

1330  NEXT  J 

1340  UN0RM=SQR<UH0RM):  T 
EST=0 

1350  FOR  1=1  TO  NPLV 
1360  X<I)=X<I>/UN0RM 
1370  TEST=TEST+X< I >*W*SG 
N<H< I >  > 

1330  NEXT  I 
1390  UNQRM=0 
1400  FOR  1=1  TO  NPLV 
1410  Z< I )=X< I )-TEST*W*SG 
N<H< I > ) 

1 420  UN0RM=UN0RM+2 < I ) *Z < 
I> 

1430  NEXT  I 

1440  IF  UNORIK IE-6  THEN 
F*=MFAIL":  RETURN  ELSE 
F$="  " 

1450  UH0RM=SQR  <  UNORM ) 
1460  FOR  I=t  TO  NPLV 
1470  2<I)=2<I)/UNORM 
1480  NEXT  I 
1490  GOSUB  740 
1500  RETURN 

1510  ’**  STRAIN  RATIO  ** 
1520  FOR  N=1  TO  NL 
1530  SUAL=SREF*SREF' < 1-E 
6  V£11AX*<E<N»  1  )*£<N>  !>•*•£ 
<N,2>*E<N,2>+E<N,3>*E<N, 
3>/2) 


1170-1260  Get  gradient  of  each 
active  constraint 


1270-1340  Sum  gradients  and 
normal i le  result 

1350-1490  Project  onto 
constant  thickness  plane.  Test 
tor  minimum  and  normalize  tinal 
resu  1  t 


1520-1560  Find  distance  trom 
tarthest  constraint  to  origin 


1540  SUAL*SQR<SUAL  > 

1550  IF  SUAL >S  THEN  S=SU 
AL 

1560  NEXT  N 
1570  RETURN 
1580  '**  IFP  ** 

1590  W=1/SQR<NPLV> 

1600  FOR  1=1  TO  NPLV 
1610  Z<I)=Ul:  H< I )=W 
1620  NEXT  I 
1630  GOSUB  7 40 
1640  S=0:  SkEF=l 
1650  GOSUB  580:  GOSUB  15 
10 

1660  FOR  1=1  TO  NPLV 
1670  H<I)=H<I)*S 
1630  NEXT  I 
1690  S=0 

1700  GOSUB  740:  GOSUB  58 

0:  GOSUB  290 

1710  RETURN 

1720  ’**  TRANSFORM  ** 

1 730  J=QXX+QVV+2*GXV : K=Q 
SS-QXV 

1740  FOR  1=1  TO  NPLV 
1750  C2=C0S<  T  < I ) ) : C2=C2* 
C2 

1760  S2=SIN<T< I) > : S2=S2* 
S2 

1770  Q 1 1 < I ) =C2*C2*GXX+S2 
*S2*QV V+2*S2*C2*  <  QXV+2*Q 
SS) 

1730  Q22< I )=S2*S2*QXX+C2 
*C2*QVV+2*C2*S2*<GXV+2*Q 
3S) 

1790  Q12(1)~(J-<Q11<1)+Q 
22<I))>/2 

1800  Q66<n  =  <J+2*K-<Qll< 

I)+Q22<I)>)/'2 

1810  NEXT  I 

1820  RETURN 

1830  ’***  INPUT  *** 

1840  CL3 

1850  PRINT  "PRESS  AHV  KE 
V  WHEN  ", "DESIRED  MATER  I 
AL", "APPEARS" 

1860  FOR  K=1  TO  750: NEXT 
1870  FOR  M=0  TO  6 
1330  IF  M=6  THEN  M$="NEU 
MATERIAL"  ELSE  GETXM, EX 
, EV , UX , ES , TPLV , XT , VT , XC , 
VC,SS,M* 

1890  CLS: PRINT  Mf: SOUND 

20,1 

1900  FOR  J=1  TO  200 
1910  IF  INKEV$<>“"  THEN 
1940 

1920  NEXT  J.M 

1930  GOTO  1850 

1940  IF  M=6  THEN  GOSUB  3 

750: GOTO  1850 

1950  CLS: PRINT  "%  ";M*J" 

T 

1960  PRINT  "HOW  MANV" 
1970  INPUT  " PLV  GROUPS"! 
NPLV 

1980  CLS  s  PR I NT " ENT  ER  PLV 
GROUP" 


1730-1810  Transformation  of 
elasticity  matrix,  assuming 
ortho tropic  laminate 
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2470  CLS: PRINT "Press  V  i 
f  printout", "of  displace 
d  result  is  desired.  Pre 
ss  N" » “if  not"; 

2480  ITER=1 

2490  FOR  1=1  TO  600: NEXT 
I 

2500  A$=INKEV$: IF  h*<>" " 
THEN  2500 

2510  CLS: RESTORE  4630 
2520  J=0:A*=INKEV$ 

2530  FOR  1=1  TO  3 
2540  READ  A*: CLS: PRINT: P 
RINT  A#: SOUND  20,1 
2550  A*=INKEV*:IF  A*="" 
THEN  2550 

2560  PRINT  Af;:FOR  KK=1 

TO  75: NEXT  KK 

2570  IF  A$=“V"  THEN  J=J+ 

l:CX<J,l>=I 

2530  NEXT  I 

2590  FOR  K=1  TO  J 

2600  ON  C:-i<K,l)  GO  SUB  32 

10,3350, 2650 , 2730 , 2390 , 2 

980,3090 

2610  L PR I NT 

2620  NEXT  K 

2630  CLS : PR I NT "F I N I SHED " 
,K* 

2640  IF  INKEV$=" "THEN264 
0  ELSE  RUN 

2650  ’  **  PLV  RATIO** 

2660  CLS: LPRINT  "Total  t 
hickness=" 

2670  LPRINT  USING  ".#### 

AAAA  ^  I*  •  JEST 

2630  LPRINT  USING  "####. 
##  P 1  i  es " ;  TEST --  TPL V 
2690  LPRINT 

2700  LPRINT  "ANGLE  RATI 
0  #F'LIES" 

2710  FOR  1=1  TO  NPLV 
2720  A=CINT< <FNDEG<T<I >> 
*1E2))/1E2 

2730  B=CINT <  <H< I )/TEST*l 
E4)  V1E4 

2740  C=CINT <  <H< I )/TPLV*l 
E2)>'1E2 

2750  LPRINT  A?TAB<6>;B;T 

AB<13>;C 

2760  NEXT  I 

2770  RETURN 

2730  ’**  STRENGTH** 

2790  LPRINT  "STRENGTH  RA 
TIOS" 

2800  LPRINT  "1=ULTIMATE 
STRAIN": 

2810  LPRINT  ">1  IS  SAFE" 
2320  FOR  1=1  TO  NL 
2830  LPRINT  "LOADING  "I 
2840  A=E<I, 1)*E<I, 1>+Evl 
,2>*E<I,2>  E<I,3>*E<I,3) 
/2 

2850  A=SCR<EMAX/A> 

2860  LPRINT  "R=”!A 
2370  NEXT 
2830  RETURN 


251 0-2620  Branch  for  variou- 
output  routines 


2650-3520  Output  routines  ana 
laminate  analysis 


2890  '**STRAINS** 

2900  LPRINT  TAB<4) * "LANI 
NATE  STRAINS" 

2910  FOR  N=1  TO  NL 
2920  LPRINT  "LOADING  "N 
2930  LPRINT  USING  "el=+# 
. ###E-03".E(N. 1>*1E3 
2940  LPRINT  USING  "e2=+# 
. ###E-03"?E(N.2>*1E3 
2950  LPRINT  USING  "*6=+# 
.  ###E-03";E<N,3>*1E3 
2960  NEXT  N 
2970  RETURN 
2980  ’  **A  MATRIX** 

2990  CLS 

3O00  LPRINT “Norm.  |A|  in 
GPa.  " 

3010  Da.l>=All:D<1.2>=A 
12: D  <  2 . 2)=A22: D<3> 3)=A66 
3020  D(1?3)=0:D(3.  1>=0:D 
<2.3>=0:D<3. 2}=0:D(2. 1 >= 
D<1. 2) 

3030  FOR  1=1  TO  3 
3040  FOR  ,T=1  TO  3 
3050  DC i , J>=DC I.  J).-'TEST/ 
1E9 

3060  NEXT  J. I 
3070  GOSUB  3430 
3080  RETURN 
3090  'A  INUERSE 
3100  LPRINT "Compliance 
(normal ized) " 

3110  LPRINT" in  1/TPa. “ 
3120  D<1> 1)=AI11:DC1.2>= 
AI12:D(2,2)=AI22:D(3,3)= 
AI66 

3130  D<1,3)=0:D(2.3>=0:D 
<3.2)=0:D<3. i>=0 
3140  D<2>  l)=Da.2> 

3150  FOR  1=1  TO  3 
3160  FOR  J=1  TO  3 
3170  D(I, T >=D ( I ,  J ) *TEST* 
1E12 

3180  NEXT  J. I 
3190  GOSUB  3430 
3200  RETURN 

3210  LPRINT  "Material  Pr 
operties" 

3220  GETXM,  EX. EV.  UX.- ES,  T 
PLV. XT . VT . XC , VC. SS . M$ 
3230  LPRINT  M$ 

3240  LPRINT  "EX=" 5  EX; "GP 
a" 

3250  LPRINT  "EV="5EV? "GP 
a" 

3260  LPRINT  "ES=";ES;"6P 

3  " 

3270  LPRINT  "UX="!UX 
3230  LPRINT  "X=" 5 XT; "MPa 

M 

3290  LPRINT  "X’="?XC;"MP 
a" 

3300  LPRINT  "V="; VT; "MPa 

M 

3310  LPRINT  "V’=" : VC. "MP 
a" 

3320  LPRINT  "S=";SS;"MPa 

M 

3330  LPRINT  "Ply  Thickne 
ss"!TPLV"m" 

3340  RETURN 


0.350  ’  **LOADS** 

3360  FOR  1=1  TO  NL 
3370  L PRINT  "LOADING  "5 1 
3330  FOR  J=1  TO  3 
3390  A$ =STR$  <  J >  s IF J =3THE 
N  A$=“  6" 

3395  A*CINT(XN< I >  J>/1E3) 
/1E3 

3400  L PRINT  "N";A*5 "="5A 
?U$ 

3410  NEXT  J,I 
3420  RETURN 
3430  ’ FANCV 

3440  LPRINT  "  t - 1 — 

- 1 - , " 

3450  LPRINT  USING  "|###. 
###" ’ D< 1 » 1 >  >  DC  1 >  2>  >  DC  1 >  3 

■) 

3460  A*="  I - t- - 

- ,  ■■ 

3470  LPRINT  A* 

3480  LPRINT  USING  "!###. 
#44" ; D<2> 1 ) , DC 2- 2 ) j  D  C  2  >  3 
) 

3490  LPRINT  h* 

3580  LPRINT  USING  "|###. 
###" J  DC3. 1 >  >  D<3>  2) i D<3>  3 
■> 

3510  LPRINT  "» - * — 

_ i _ i 

3520  RETURN 


3750  PR I NT "REVIEW  OR  NEW 

3760  INPUT  "DATA  (p/NV; 
A$ 

3770  IF  A*="R"  then  3950 
3780  PRINT"WHICH  MATERIA 
L"! PRINT-WILL  VOU" 

3790  INPUT  "REPLACE  <0-5 
)  "5 1 

3800  INPUT  "EX<GPa’>=":EX 
3S10  INPUT  "EV<GPa)="5EV 
3820  INPUT  "UX=";UX 
3330  INPUT  "ESCGPa)=";ES 
3340  INPUT  "XCHPa)=M;X 
3S50  INPUT  "X’ <MPa)*"«XX 
3860  INPUT  "VCriPa>  =  "; V 
3870  INPUT  "V’ <MPa>="; VV 
3830  INPUT  "S<MPa)=";S 
3390  If  (PUT  "THICKNESS  Cm 
)  =  " >  TPLV 

3900  INPUT  "NAME  <15  CHR 

MAV  S  H  !  Mf 

3910  PUT5iI,EX,EV,UX,ES,T 
PLV»X» V>XX» VV>S>M$ 

3920  PRINT  "ADDITIONAL": 
INPUT  "CHANGES  <V/N>" ’ AT 
3930  IF  A*="V"  THEN  375A 
3940  RETURN 

3950  PR I NT" REVIEW  WHICH" 

: INPUT "MATERIAL  <0-5> " : m 
3960  G0SUB32 1 0 
3970  GOTO  3920 


3980  OPEN  "I" >41, “CASltD 
AT  A" 

3990  FOR  1=0  TO  5 
4000  INPUT  #1,EX,EV,UX,E 
S»T»X> V>XX>VV>S>M$ 

4010  PUTXI >  EX>  EV, UX>  E3>  T 
,X,V,XX,  W,S,M* 

4020  NEXT 
4030  CLOSE  #1 
4040  DELETE  45 
4060  ’  ***RGTAT  I  ON*'** 

4065  CLS: PRINT  “WORK I NO" 
4070  A=-ROT-ROTSUM: ROTSU 
M=ROTSUM+A 
4080  A=FNRADOA) 

4090  C=C0SCA>:C2=C*C 
4100  S=S INCA)!  S2=S*S 
4110  FOR  1=1  TO  NL 
4120  R(1 )=XN( I , 1 >*C2+XN( 
I . 2>*32+XN< I , 3>*2*S*C 
4130  R  < 2 )=XN  < I , U *S2+XN C 
I , 2)*C2-XN( I , 3>*2*S*C 
4140  R  < 3 >=-XN <  I  >  1  > *C*S+X 
N <  1 , 2 ) *C*S+XN<  I>  3>*C C2-S 
2) 

4150  XN< I > 1 >=R< 1 > : XN< I , 2 
>=R<2>  SXNC I >  30=R<3) 

4160  NEXT  I 
4170  RETURN 

4180  ’ ** ANGLE  SEARCH** 
4190  INFUT  "LOWER  SEARCH 
LIMIT" “SI 

4195  INFUT  "UFPER  SEARCH 
LIMIT" ?  S2 

42O0  IF  S2<=S1  THEN  4190 
4210  INPUT  " MAX.  ERROR"? 
El 

4220  ROTSUM=0 
4230  U< 1 >=S1 : U(5)=S2 
4240  U(3>=<Sl+S2>/2 
4250  U<2>  =  <U<3>HK1)V2 
4260  U<4)*<U<5)+U<3) )/'2 
4265  AA=1?BB=5 
4270  FOR  I I=AA  TO  BB 
4280  ROT=U< 1 1 ) 

4290  QOSUB  4O60?  GOSUB  46 
50 

4300  U < 1 1 >  =H 

4310  NEXT 

4320  UMIN=U<l)iJ=l 

4330  FOR  1=2  TO  5 

4340  IF  U< I ) CUMIN  THEN  U 

MIN=U( I >  s  J=I 

4350  NEXT 

4360  IF  J=1  THEN  4820 
4365  IF  J=5  THEN  4770 
4370  V<1)=U< J-1>:XC1>=UC 
J-U 

4388  VC 2) =U < J > ! X ( 2 >  =U< J > 
4390  V<3)=UC J+l ) sXC3>=UC 

J+l) 

4400  U<1>=X<  1  >!U(3;,=X<2) 
:UC5)=X<3) 

4410  U<1>=V<1 j!UC3>=VC2) 
!'K5)*V<3> 

4420  IF  U<5)-U< 1 ><=E1  TH 

EN  4540 


3y80-4040  Routine  For 
automatically  entering  material 
properties  From  cassette  tape. 
To  use,  line  45  should  read 
GOTO  3980,  and  run  program  with 
tape  still  connected 


4070-4170  TransForm  loads  to 
new  axis  system. 


4220-4610  One-dimensional 
search  For  best  rigid  body 
rotation 


m 


1990  PRINT  "ORIENTATIONS 

2000  FOR  I*i  TO  350 
2010  NEXT 
2020  CLS 

2030  FOR  1=1  TO  NPLV 
2040  PRINT  "PLV  "51 
2050  INPUT  T<I> 

2060  T<I)=FNRAD<T<I>) 
2070  NEXT  I 

2030  INPUT  "NUMBER  OF  LO 
ADINGS=";NL 
2090  FOR  1=1  TO  NL 
2100  CLS: PRINT  "LOADING 
I 

2110  INPUT  "Nl(MPa>=";XN 
<1,1) 

2120  INPUT  "N2(MPa)=";XN 
CI»2) 

2130  INPUT  "N6<MPa)=";XN 

<  1, 3) 

2140  FOR  J=1  TO  3 

2150  XN< I , J)=XN< I , J)*1E6 

2160  NEXT  J: NEXT  I 

2170  RETURN 

2130  ’*■*  INUARIENTS*'* 

2190  UV=l/< 1 -UX*UX+EV/EX 
) 

2200  QXX=UV*EX* 1 E9 : QVV=U 
Y*EV*1E9 

2210  QX V=UV*UX*EV* 1 E9 : QS 

S=ES*1E9 

2220  U<1)=XT/EX 

2230  U<2)=XC/EX 

2240  U<3)»VT/'EV 

2250  U<4)=VC/EV 

2260  U  <  5 )  =SS  •"ES-'SQR  <  2 ) 

2270  EMAX=U  <  1  )*U<  1 )  ✓  1 E6 

2230  FOR  1=2  TO  5 

2290  U< I )=U< I >*U< I )/lE6 

2300  IF  U< I XEMAX  THEN  E 

MAX=U< I ) 

2310  NEXT  I 
2320  RETURN 
2330  ’**  OUTPUT  ** 

2340  SOUND  15,2:SOUND50, 

2 

2350  K$="Hit  any  key  to 
cont.  ":U$="MN/'m" 

2360  CLS:TEST=0 
2365  IF  F#="ROT"  THEN  LP 
RINT  "RIGID  BODY” , "ROTAT 
ION  OF" ;U<3>, "DEGREES" 
2370  FOR  1=1  TO  NPLV 
2380  TEST=TEST+H< I > : NEXT 
I 

2390  PRINT  "TOTAL  THICKN 
ESS=" 

2400  PRINT  TEST! "  m.  " 

2410  PRINT  USING  "####.# 

#  Plies" ; TEST^TPLV 
2420  PRINT  K*; 

2430  LOCATE  0,0 
2440  A#=INKEV$: IF  A$<>"" 
THEN  2440 

2450  IF  INKEV*=H"  THEN  2 
450 


2190-2210  Calculate  Q's  -from 
engineering  constants 


2220-2310  Find  smallest  strain 
component.  It  becomes  the 
radius  o-f  the  strain-sphere 


4430  U<2>  =  <lK3>+U<n>/2 
4440  LK  4  > = <  IK  5  >  +  LK  3 } > / 2 
4450  FOR  11=2  TO  4  STEP 
2 

4460  ROT=U<  1 1  !> 

4470  GOSUB  4060: GOSUB  46 
50 

4430  U(II)=H 
44?0  NEXT 
4500  GOTO  4320 
4540  R0T=LK3> 

4550  GOSUB  4060: GOSUB  46 
50 

4560  ROT=0: GOSUB  4060 
4570  PRINT  "OPT.  ORIENT A 
TI0H="?U<3) 

4580  PRINT  "TOTAL  OF  " ? U 
<3)/TPLVj "PLIES"! 

4590  SOUND  1 5>  2: SOUND  30 
,  2 

4600  IF  INKEV$=" "  THEN  4 
600 

4605  Ft* "ROT" 

4610  GOTO  2330 
4620  DATA  5E-5,.  1, IE-6 
4630  DATA  Ply  properties 
, Loads  >  Tota 1  th  i ck ness  fe 
ply  r at i os » Strength 
ratios 

4640  DATA  Laminate  strai 
ns>  Sti f fness  matr ix> Comp 
lienee  matrix 

4650  ’ **0PT.  RATIO** 

4660  GOSUB  1530 

4670  GOSUB  1160 

4680  IF  FS= " FAIL"  THEN  4 

723 

4690  GOSUB  340 

4700  IF  F$=“FAIL"  THEN  4 

720 

4710  GOTO  4673 
4720  H=0 

4730  FOR  1=1  TO  NPLV 
4740  H=H+H(I) 

4750  NEXT  I 
4760  RETURN 

4770  U  <  1 )  =IJ<  4  > :  U  <  1 )  =U  4  ) 
: U  <  2  > =U<  5 ) : U  < 2 ) =U  <  5  > 

4730  DEL=IJ(2>-U<  1  > 

4790  U < 3 )  =U ( 2) +DEL :  LK 4 )  = 

U  <  3  )  +DEL  :IJ<5  >  =U  <  4  )  +DEL 

48O0  AA=3: BB=5 

4810  GOTO  4270 

4320  U  <  5  >  =U  <  2 ) : U  < 5 ) =U  C 2 ) 

:U<4)*U< 1 ) : U<4)=U< 1 > 

4330  DEL=IJ  <  5  >  -U  <  4 ) 

4840  U<3)=LK4)-DEL:U<2)= 
IJ  3  >  -DEL :  U  <  1 )  =U  <  2 i  -DEL 
4350  AA=1 : BB=3 
4860  GOTO  4270 


4650-47*0  Vers i  on  of  MAIN  used 
by  op  t  i  murri  angle  search  to 
optimize  ply  ratios  at  each 
trial  or i en  tat i on 


4770—4860  Routine  to  adjust 
search  bounds  it  minimum  not 
found  in  given  limits 
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